Clutter suppression in ultrasonic imaging systems

ABSTRACT

A method of ultrasound imaging is provided herein. The method includes the following stages: transmitting ultrasound radiation towards a target and receiving reflections of the ultrasound radiation from a region of the target in a main reflected signal and one or more auxiliary reflected signals, wherein each one of the reflected signals comprises an input dataset and is associated with a beam with a different and distinct beam pattern; compounding the input datasets from the main reflected signal and one or more auxiliary reflected signals, by the use of a compounding function, said compounding function using parameters derived from spatial analysis of the input datasets.

This application is a continuation of Ser. No. 13/408,136 filed Feb. 29,2012, which is a continuation-in-part of U.S. patent application Ser.No. 13/234,927 filed Sep. 16, 2011, which is a continuation-in-part ofU.S. patent application Ser. No. 11/722,246 filed on Jun. 20, 2007, nowU.S. Pat. No. 8,045,777, issued Oct. 25, 2011, which is a national stageentry of PCT/IL2005/001383 filed on Dec. 27, 2005 which claims priorityto U.S. Provisional Application Ser. No. 60/640,368 filed Dec. 30, 2004,all of which are incorporated herein by reference in their entirety.

FIELD OF THE INVENTION

The present invention relates generally to ultrasonic imaging systems,e.g., for medical imaging, and particularly to methods and systems forsuppressing clutter effects in ultrasonic imaging systems.

BACKGROUND OF THE INVENTION

Ultrasonic medical imaging plays a crucial role in modern medicine,gradually becoming more and more important as new developments enter themarket. One of the most common ultrasound imaging applications isechocardiography, or ultrasonic imaging of the cardiac system. Otherwidespread applications are obstetrics and gynecology, as well asabdominal imaging, to name a few. Ultrasonic imaging is also used invarious other industries, e.g., for flaw detection during hardwaremanufacturing.

Ultrasonic imaging systems typically produce relatively noisy images,making the analysis and/or diagnosis of these images a task for highlytrained experts. One of the most problematic imaging artifacts isclutter, i.e., undesired information that appears in the imaging plane,obstructing data of interest.

One of the main origins of clutter in ultrasonic imaging is effectiveimaging of objects outside the probe's mainlobe, also referred to assidelobe clutter. For example, in echocardiography, the dominantreflectors outside the probe's mainlobe are typically the ribcage andthe lungs.

Another origin of clutter is multi-path reflections, also calledreverberations. In some cases, the geometry of the scanned tissue withrespect to the probe, as well as the local reflective characteristics ofthe tissue, causes a noticeable percentage of the transmitted energy tobounce back and forth in the tissue before reaching the probe. As aresult, the signal measured for a specific range with respect to theprobe may include contributions from other ranges, in addition to thedesired range. If the signal emanating from other ranges is caused byhighly reflective elements, it may have a significant effect on theimage quality.

A common method for enhancing the visibility of the desired ultrasonicimage relative to the clutter, particularly in patients with lowechogenicity (a common phenomenon among obese patients), isadministering contrast agents. Such agents enhance the ultrasonicbackscatter from blood and aid in its differentiation from thesurrounding tissue. This method is described, for example, by Krishna etal., in a paper entitled “Sub-harmonic Generation from UltrasonicContrast Agents,” Physics in Medicine and Biology, vol. 44, 1999, pages681-694, which is incorporated herein by reference.

Using harmonic imaging instead of fundamental imaging, i.e.,transmitting ultrasonic signals at a certain frequency and receiving attwice the transmitted frequency, also reduces clutter effects. Spenceret al. describe this method in a paper entitled “Use of Harmonic Imagingwithout Echocardiographic Contrast to Improve Two-Dimensional ImageQuality,” American Journal of Cardiology, vol. 82, 1998, pages 794-799,which is incorporated herein by reference.

U.S. Pat. No. 6,251,074, by Averkiou et al., issued on Jun. 26, 2001,titled “Ultrasonic Tissue Harmonic Imaging,” describes an ultrasonicdiagnostic imaging system and methods, which produce tissue harmonicultrasonic images from harmonic echo components of a transmittedfundamental frequency. Fundamental frequency waves are transmitted by anarray transducer to focus at a focal depth. As the transmitted wavespenetrated the body, the harmonic effect develops as the wave componentsbegin to focus. The harmonic response from the tissue is detected anddisplayed, while clutter from fundamental response is reduced byexcluding fundamental frequencies.

Furthermore, image-processing methods have been developed for detectingclutter-affected pixels in echocardiographic images by means ofpost-processing. Zwirn and Akselrod present such a method in a paperentitled “Stationary Clutter Rejection in Echocardiography,” Ultrasoundin Medicine and Biology, vol. 32, 2006, pages 43-52, which isincorporated herein by reference.

Other methods utilize auxiliary receive ultrasound beams. In U.S. Pat.No. 8,045,777, issued on 25 Oct. 2011, titled “Clutter Suppression inUltrasonic Imaging Systems,” Zwirn describes a method for ultrasonicimaging, comprising: transmitting an ultrasonic radiation towards atarget; receiving reflections of the ultrasonic radiation from a regionof the target in a main reflected signal and one or more auxiliaryreflected signals, wherein each one of the reflected signals isassociated with a different and distinct beam pattern, wherein all ofthe reflected signals have an identical frequency; determining ade-correlation time of at least one of: the main reflected signal andthe one or more auxiliary reflected signals; applying a linearcombination to the main reflected signal and the one or more auxiliaryreflected signals, to yield an output signal with reduced clutter,wherein the linear combination comprises a plurality of complex numberweights that are being determined for each angle and for each rangewithin the target tissue, wherein each complex number weight is selectedsuch that each estimated reflection due to the clutter is nullified,wherein a reflection is determined as associated with clutter if thedetermined de-correlation time is above a specified threshold.

U.S. patent application 2009/0141957, by Yen and Seo, published on Jun.4, 2009, titled “Sidelobe Suppression in Ultrasound Imaging using DualApodization with Cross-Correlation,” describes a method of suppressingsidelobes in an ultrasound image, the method comprising: transmitting afocused ultrasound beam through a sub-aperture into a target andcollecting resulting echoes; in receive, using a first apodizationfunction to create a first dataset; in receive, using a secondapodization function to create a second dataset; combining the twodatasets to create combined RF data; calculating a normalizedcross-correlation for each pixel; performing a thresholding operation oneach correlation value; and multiplying the resulting cross-correlationmatrix by the combined RF data.

An additional class of currently available methods for handling clutteris a family of clutter rejection algorithms, used in color-Doppler flowimaging. These methods estimate the flow velocity inside cardiacchambers or other blood vessels and suppress the effect of slow-movingobjects; assuming that the blood flow velocity is significantly higherthan the motion velocity of the surrounding tissue. These methods aredescribed, for example, by Herment et al. in a paper entitled “ImprovedEstimation of Low Velocities in Color Doppler Imaging by Adapting theMean Frequency Estimator to the Clutter Rejection Filter,” IEEETransactions on Biomedical Engineering, vol. 43, 1996, pages 919-927,which is incorporated herein by reference.

SUMMARY OF THE INVENTION

Embodiments of the present invention provide methods and devices forreducing clutter effects in ultrasonic imaging systems.

The present invention provides a method of ultrasound imaging, saidmethod comprising: transmitting ultrasound radiation towards a targetand receiving reflections of the ultrasound radiation from a region ofthe target in a main reflected signal and one or more auxiliaryreflected signals, wherein each one of the reflected signals comprisesan input dataset and is associated with a beam with a different anddistinct beam pattern; compounding the input datasets from the mainreflected signal and one or more auxiliary reflected signals, by the useof a compounding function, said compounding function using parametersderived from spatial analysis of the input datasets, said compoundingfunction parameters being derived from at least one of:

(a) The local phase difference and/or magnitude difference and/ormagnitude ratio and/or complex signal ratio between different receivebeams, and/or functions of one or more of the aforementioned parameters;(b) Spatial functions of the local phase difference and/or magnitudedifference and/or magnitude ratio and/or complex signal ratio betweendifferent beams;(c) The local difference and/or ratio between spatial functions of thelocal magnitude and/or phase and/or complex signal in the differentreceive beams;(d) Applying multiple temporal band-pass filters to the local magnitudeand/or phase and/or complex signal in the different receive beams,applying spatial analysis to the outputs of multiple temporal band-passfilters and compounding the results;(e) Applying multiple temporal band-pass filters to the local phasedifference and/or magnitude difference and/or magnitude ratio and/orcomplex signal ratio between different receive beams, applying spatialanalysis to the outputs of multiple temporal band-pass filters andcompounding the results; or(f) Spatial-temporal functions of the local phase difference and/ormagnitude difference and/or magnitude ratio between different beams,wherein spatial derivatives of the local phase difference and/ormagnitude difference and/or magnitude ratio between different beams, areapplied after local temporal filtering; or wherein the local differenceand/or ratio between spatial functions of the magnitude and/or phaseand/or complex signal in the different beams are applied after localtemporal filtering.Other aspects of the present invention are detailed in the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention for clutter suppression in ultrasonic imaging systems isherein described, by way of example only, with reference to theaccompanying drawings.

With specific reference now to the drawings in detail, it is stressedthat the particulars shown are by way of example and for purposes ofillustrative discussion of the preferred embodiments of the presentinvention only, and are presented in the cause of providing what isbelieved to be the most useful and readily understood description of theprinciples and conceptual aspects of the invention. In this regard, noattempt is made to show structural details of the invention in moredetail than is necessary for a fundamental understanding of theinvention, the description taken with the drawings making apparent tothose skilled in the art how the several forms of the invention may beembodied in practice.

FIG. 1A is a schematic, pictorial illustration of an ultrasonic imagingsystem, in accordance with an embodiment of the present invention;

FIG. 1B is a schematic, pictorial illustration of a probe used in anultrasonic imaging system, in accordance with an embodiment of thepresent invention;

FIG. 2 is a schematic illustration of the beam pattern of two exemplaryreceive beams, defined in terms of gain as a function of the azimuthangle, in accordance with an embodiment of the present invention;

FIG. 3A is a schematic illustration of the expected magnitude ratiobetween the two receive beams of FIG. 2, defined for a point-likereflector as a function of the azimuth angle, wherein the azimuth scaleis different than that in FIG. 2, providing a zoom on angles close to0°, in accordance with an embodiment of the present invention;

FIG. 3B is a schematic illustration of the expected phase differencebetween the two receive beams of FIG. 2, defined for a point-likereflector as a function of the azimuth angle, wherein the azimuth scalematches that of FIG. 3A, in accordance with an embodiment of the presentinvention.

DETAILED DESCRIPTION OF EMBODIMENTS System Description

In broad terms, the present invention relates to methods and systems forsuppressing clutter effects in ultrasonic imaging systems. In someembodiments, these methods and systems also result in enhancement of theultrasound images' lateral resolution, i.e., the resolution alongdirections perpendicular to the centerline of an ultrasound beam.

Before explaining at least one embodiment of the invention in detail, itis to be understood that the invention is not limited in its applicationto the details of construction and the arrangement of the components setforth in the following description or illustrated in the drawings. Theinvention is capable of other embodiments or of being practiced orcarried out in various ways. Also, it is to be understood that thephraseology and terminology employed herein is for the purpose ofdescription and should not be regarded as limiting.

FIG. 1A is a schematic, pictorial illustration of an ultrasonic imagingsystem 20, in accordance with an embodiment of the present invention.

System 20 comprises an ultrasound scanner 22 that scans organs of apatient using ultrasound radiation. A display unit 24 displays thescanned images. A probe 26, connected to scanner 22 by a cable 28, istypically held against the patient body in order to image a particularbody structure, such as the heart (referred to as a “target”).Alternatively, the probe may be adapted for insertion into the body,e.g., in transesophageal or transvaginal configurations. The probetransmits and receives ultrasound beams required for imaging. Scanner 22comprises control and processing circuits for controlling probe 26 andprocessing the signals received by the probe.

FIG. 1B is a schematic, pictorial illustration of probe 26 used inimaging system 20, in accordance with an embodiment of the presentinvention. The probe comprises an array of transducers 30, e.g.,piezoelectric transducers, which are configured to operate as a phasedarray. On transmission, the transducers convert electrical signalsproduced by scanner 22 into a beam of ultrasound radiation transmittedinto the patient's body. On reception, the transducers receive theultrasound radiation reflected from different body tissues, and convertit into electrical signals sent to scanner 22 for processing.

Probe 26 typically comprises several dozens and up to several hundredsof transducers 30 arranged in a horizontal linear array. The horizontalaperture of the array, typically of the order of several centimeters,affects the minimal horizontal beam-width of the probe and the resultinghorizontal angular resolution. The vertical beam-width may be adjustedby an acoustic lens. Some probes, referred to as “1½ dimensionalprobes,” comprise several rows of transducers in the vertical dimension,providing a vertical sector-like beam pattern. Other probes comprise acomplete two-dimensional (or multi-dimensional) array of transducers,enabling control over both horizontal and vertical directional beampatterns. The terms “horizontal” and “vertical” are used here solely forconvenience, as the array may be positioned during imaging in anyappropriate orientation relative to the patient's body.

A single array of transducers may generate different beam patterns,whose mainlobe may be pointed at different directions. On transmission,the direction of the mainlobe is typically set by adjusting the relativephase and/or the relative time delay of the signals fed to the differenttransducers. On reception, the direction of the mainlobe is usually setby adjusting the phase shift and/or time delay introduced to the signalreceived by each transducer.

The beam pattern on both transmission and reception may be adjusted byapplying different apodizations. On reception, apodization is a processof multiplying the output signal of each transducer by a multiplicativecoefficient (“apodization coefficient”), before combining the outputs ofall transducers to generate an overall array output signal. A reciprocalprocess is performed on transmission (in many cases a constantapodization is used on transmission).

Let k be the transducer index (k should go over all transducers even ifthe transducer array comprises of more than one dimension), s_(k) be thesignal measured by transducer k, a_(k) be the apodization coefficient oftransducer k on reception, and φ_(k) be the phase shift for transducer kon reception. Beamforming on reception, for obtaining the signal S attime t, may be performed using eq. (1):

S(t)=Σ_(k) a _(k)(t)e ^(jφ) ^((t)) s _(k)(t)  (1)

Alternatively, when using time-delays instead of phase shifts, whereτ_(k) is the time-delay for transducer k, one may use eq. (2) forbeamforming on receive:

S(t)=Σ_(k) a _(k)(t)s _(k)(t−τ _(k))  (2)

Similar equations may be used for combined time-delay and phase-shiftdesign. Comparable equations may also be utilized on transmission.

Depending on system design, the signal reaching the differenttransducers may be sampled (by one or more analog-to-digital convertersin scanner 22) at different beamforming stages, e.g., for eachtransducer separately, for different groups of transducers, or for alltransducers combined. Eq. (1) and eq. (2) and combinations or variationsthereof are applicable to all these cases.

Moreover, some systems perform sampling for all transducers combined(i.e., at the beam level), but use more than one set of apodizationcoefficients a_(k) and/or more than one set of phase shifts φ_(k) and/ormore than one set of time-delays τ_(k), a setting commonly referred toas multi-line acquisition, or MLA. In MLA configurations, the beampattern used on transmission is typically wider than those used onreception, so as to provide ultrasound energy to most or all of thevolume covered by the different receive beams.

Certain scanners 22 apply filtering matched to the transmitted waveformand/or down-conversion analogically, before sampling, whereas othersapply matched filtering and/or down-conversion digitally. The matchedfiltering and/or down-conversion may be applied in any of thebeamforming stages.

Some scanners 22 sample the real data at one or more of the beamformingstages, whereas others use quadrature receivers, providing complexsamples, including both a real in-phase component and an imaginaryquadrature component.

Further types of scanners 22 use multiple concurrent beams ontransmission, each of which having a different waveform (e.g., differentcentral transmission frequency, different pulse encoding configuration).In these cases, each receive beam is typically generated using aspecific matched filter, corresponding to one of the transmissionwaveforms.

Note that in some probe designs, a planar array is used instead of aphased array. In such cases, only a small set of predefined time-delaysand/or phase-shifts is supported by the transducer array, wherein thebeam formed by each such set is typically sampled separately. Furtherprobe designs do not use transducer arrays, but instead allocate asingle transducer to each transmit and/or receive beam.

Probe 26 typically scans the region of interest, which may beone-dimensional (1D), two-dimensional (2D) or three-dimensional (3D), byone or more of the following:

(a) Electronic beam-steering, i.e., changing the direction of themainlobe by adjusting the relative phase and/or the relative time delayof the signals fed to the different transducers.(b) Electronic change of the phase center of the transducer array,performed by adjusting the apodization pattern, e.g., by setting certainregions of the apodization pattern to 0. As a result, the vector betweenthe phase center of the beam and the peak of its mainlobe (the “beamdirection”) may be spatially shifted, producing another vector parallelto the former one.(c) Mechanical beam-steering, i.e., mechanically moving portions of orall of the transducer array, thus affecting the direction of itsmainlobe.

System Configurations for Clutter Suppression

In embodiments of the present invention, one or more auxiliary receivebeams are used for suppressing clutter effects. The auxiliary receivebeams may be generated using any configuration known in the art,including but not limited to those described hereinabove.

Auxiliary receive beams may be allocated separately for each beam usednormally for imaging (“main beam”), or for groups of main beams. Thelatter configuration may be useful, e.g., in MLA configurations, wheremultiple main beams are typically used concurrently on reception,corresponding to a single transmission beam, and one or more auxiliarybeams may be defined.

The auxiliary receive beams may be generated at the same time as thecorresponding main receive beam or beams, or at different times.Combinations can also be considered, i.e., some of the auxiliary receivebeams may be generated at the same time as the corresponding mainreceive beam or beams, and others may be generated at different times.Various sequences of relative timing between the different beams may bedefined, e.g., alternating generation of main and auxiliary receivebeams. If some or all of the auxiliary receive beams are generated atdifferent times than the corresponding main receive beam or beams,separate transmission beams may be used for the main and auxiliaryreceive beams or groups thereof, using the same or different beampatterns.

In exemplary embodiments, one main receive beam and one auxiliaryreceive beam are utilized, corresponding to a single transmission beamor two separate transmission beams. In such a case, the auxiliaryreceive beam may have a different beam pattern than the main receivebeam, e.g., it may be pointed at the same direction but its beam widthmay be wider (or narrower).

In further exemplary embodiments, corresponding to an MLA configuration,multiple concurrent main receive beams (e.g., 16 or 32) may be used witha single transmission beam, and one auxiliary receive beam is added, tobe used for all concurrent receive beams. In most MLA cases, the beampatterns for all main receive beams are similar, but each main receivebeam points at a slightly different direction. In order to providereference information for all main receive beams, the auxiliary receivebeam may be wider than some of them, e.g., its width may match that ofthe transmission beam.

In even further exemplary embodiments, multiple main receive beams maybe used with a single transmission beam, wherein one or more auxiliaryreceive beams may be added for each main receive beam or for each groupof main receive beams.

In certain embodiments, one or more receive beams may be treated as bothmain receive beams and auxiliary receive beams, in the context ofdifferent computations.

Clutter Suppression Method

In the embodiments of the present invention, compound information isprovided by one or more main receive beams and one or more auxiliaryreceive beams so as to provide clutter suppressed outputs. The processmay be applied in various processing phases of scanner 22, and todifferent types of data, e.g., analogically or digitally, to real orcomplex samples, before or after matched filtering and/or down-sampling.

In some cases, conversions which would be easily understood by peopleknowledgeable in the art may be required. For example, some embodimentsrequire that the input data for each receive beam would be complex. Insuch cases, if the input data for one or more receive beams is real, aHilbert transform may be applied to the real data for each receive beam,and an inverse Hilbert transform may be applied to the cluttersuppression outputs if necessary.

In configurations where scanner 22 performs sampling for each transducerseparately or for different groups of transducers, the cluttersuppression process may be applied separately for different groups oftransducers. In such cases, “synthetic” main and auxiliary receive beamsmay be defined for computation purposes only, for each group oftransducers separately, and the clutter suppressed outputs for thedifferent transducer groups may be used for further stages ofbeamforming.

In certain embodiments, the information for each main receive beamand/or auxiliary receive beam may be described as an array of real orcomplex measurements, each of which corresponding to a certain volumecovered by the mainlobe (and sidelobes) of the applicable beam, betweenconsecutive iso-time surfaces of the ultrasound wave within the medium(with respect to the probe 26), typically but not necessarily matchingconstant time intervals. Each such volume is commonly referred to as avolume pixel, or voxel. The samples are commonly referred to asrange-gates, since the speed of sound does not change significantlywhile traversing soft tissues, so that iso-time surfaces canapproximately be referred to as iso-range surfaces.

In some embodiments, a 1D, 2D or 3D region of interest is scanned byscanner 22. The scanning may be performed by any method known in theart, including but not limited to the techniques mentioned hereinabove.For example, different beams may have the same phase center butdifferent beam directions. In such cases, a polar coordinate system istypically used in 2D scanning, and a spherical coordinate system istypically used in 3D scanning. When using a polar coordinate system, thelocation of each voxel may be defined by the corresponding range-gateindex and the angular direction of the beam with respect to thebroadside of probe 26, wherein the probe's broadside is defined by aline perpendicular to the surface of probe 26 at its phase center, andwherein said angular direction may be defined in a Euclidian space by anazimuth angle, or by the u coordinate in sine-space. Similarly, whenusing a spherical coordinate system, the location of each voxel may bedefined by the corresponding range-gate index and the angular directionof the beam with respect to the broadside of probe 26, wherein the angledirection may be defined either by the azimuth and elevation anglesand/or by the (u,v) coordinates in sine-space.

Additionally or alternatively, different beams may have different phasecenters, wherein the beams may or may not be parallel to each other. Inthe case of parallel beams, a 2D or 3D Cartesian coordinate system maybe used, wherein the location of each voxel is defined by itscorresponding range-gate index and the location of the phase center.

Regardless of the scanning method and coordinate system used, thedataset collected by one or more receive beams, which may or may not beconcurrent, in each time-swath (“frame”) may be organized in a 1D, 2D or3D array (“scanned data array”), using any voxel arrangement known inthe art, wherein each index into the arrays relates to a different axis(e.g., in a polar coordinate system, a range-gate index and an azimuthindex may be utilized), so that voxels which are adjacent to each otherin one or more axes of the coordinate system also have similar indicesin the corresponding axes. In certain cases, the scanned data array mayuse a coordinate system which does not correspond to the scanningmethod, e.g., all receive beams may have the same phase center butdifferent beam directions, but the scanned data array may utilizeCartesian rather than polar or spherical coordinates. This may beobtained by spatial interpolation and/or extrapolation, using any methodknown in the art, e.g., nearest neighbor interpolation, linearinterpolation, spline or smoothing spline interpolation, and so forth.

In some embodiments, the scanned data array adheres to the followingrequirements, in which case the scanned data arrays are referred to as“common configuration scanned data arrays”:

(a) All the data in a scanned data array corresponds to receive beamswith the same beam pattern, except, perhaps, for translations and/orrotations.(b) All the data in a scanned data array corresponds to receive beamsusing the same central frequency on reception.(c) All the data in the scanned data array corresponds to receive beamsused in conjunction with transmit beams using the same waveforms and thesame transmit beam pattern, except, perhaps, for translations and/orrotations.(d) The data in a scanned data array originates from receive beams whichmay be concurrent, divided into multiple groups each of which isconcurrent (this is the most common case when using MLA configurations),or non-concurrent. In other embodiments, one or more of requirements(a)-(d) may not be met.

Compounding Functions

In exemplary embodiments, the compounding of information provided by amain receive beam and one or more associated auxiliary receive beams isperformed by computing a function of the local signal in two or morereceive beams (“compounding function”), wherein the term local signalrefers to a certain volume region, which is partially or wholly coveredby more than one of the receive beams, e.g., a certain range-gate indexof the different receive beams in cases where multiple receive beamshave the same sampling frequency. The local signal in each such receivebeam may be generated by setting the matched filters to fundamentalimaging, i.e., using a reception frequency band centered at the centerfrequency of the corresponding transmission beam, or to harmonicimaging, i.e., using a reception frequency band centered at a frequencywhich is a natural number multiplied by the center frequency of thecorresponding transmission beam (harmonies of the transmissionfrequencies). In some embodiments, different receive beams may usedifferent harmonies of the transmission frequencies, including but notlimited to first harmonies (i.e., fundamental imaging). The output ofthe compounding function for each entry into the scanned data arrayand/or each frame may be used instead of the measured signal for themain receive beam in any further processing performed by ultrasoundimaging systems, such as:

(a) Log-compression, i.e., computing the logarithm of the local signalmagnitude.(b) Time-gain control, i.e., applying a range-dependent correction forthe signal attenuation within the medium.(c) Global gain control.(d) Doppler processing, e.g., for color-Doppler flow imaging,tissue-Doppler imaging, and/or pulse-Doppler studies. Note that Dopplerprocessing typically entails analyzing data from two or more pulses foreach applicable voxel.

In some embodiments, the compounding function may use parameters(“compounding function parameters”) derived from at least one of:

(a) The local phase difference and/or magnitude difference and/ormagnitude ratio between different receive beams, and/or functions of oneor more of the aforementioned parameters, e.g., the magnitude ratioconverted into decibel units.(b) Spatial functions of the local phase difference and/or magnitudedifference and/or magnitude ratio between different beams, e.g., variousspatial derivatives of said parameters, using any operator known in theart, such as

$\frac{1}{3}\begin{pmatrix}{- 1} & 0 & 1 \\{- 1} & 0 & 1 \\{- 1} & 0 & 1\end{pmatrix}\mspace{14mu} {and}\mspace{14mu} \frac{1}{3}{\begin{pmatrix}{- 1} & {- 1} & {- 1} \\0 & 0 & 0 \\1 & 1 & 1\end{pmatrix}.}$

(c) The local difference and/or ratio between spatial functions of thelocal magnitude and/or phase and/or complex signal in the differentreceive beams, e.g., the ratio of local magnitude derivatives along acertain spatial axis; the ratio of local normalized magnitudederivatives along a certain spatial axis, wherein a local normalizedderivative is defined as the local signal derivative divided by thelocal signal; the ratio of local magnitude spatial derivatives, whereinthe spatial derivative is computed using a “Laplacian filter”,

${i.e.},{{\frac{1}{8}\begin{pmatrix}{- 1} & {- 1} & {- 1} \\{- 1} & 8 & {- 1} \\{- 1} & {- 1} & {- 1}\end{pmatrix}};}$

the ratio of local magnitude normalized spatial derivatives, wherein thespatial derivative is computed using a “Laplacian filter”; the ratio ofnormalized magnitude spans, wherein a normalized signal span is definedas the difference between the maximal and minimal value, divided by themean value, wherein an outlier rejection scheme may be applied prior tothe statistical functions minimum, maximum and/or mean; and so forth.(d) Temporal functions of the local phase difference and/or magnitudedifference and/or magnitude ratio between different beams, e.g., varioustemporal low-pass, band-pass or high-pass filters.(e) Spatial-temporal functions of the local phase difference and/ormagnitude difference and/or magnitude ratio between different beams,e.g., spatial derivatives of the local phase difference and/or magnitudedifference and/or magnitude ratio between different beams, applied afterlocal temporal filtering; local difference and/or ratio between spatialfunctions of the magnitude and/or phase and/or complex signal in thedifferent beams, applied after local temporal filtering; and so forth.

In that context, spatial filtering may be applied to the scanned dataarray for a specific frame, whereas the temporal filtering may beapplied to a specific voxel or group of adjacent voxels over multipleconsecutive frames.

In some embodiments, one or more of the following assumptions may beutilized when selecting and/or making use of the compounding functionparameters:

(a) The local phase difference and/or magnitude difference and/ormagnitude ratio are indicative of the angular direction of the signalsource. By definition, non-clutter signals originate from angulardirections close to that of the mainlobe of a main receive beam. Forexample, when using a main receive beam and an auxiliary receive, whosemainlobes point at the same spatial angle, one would expect that forreflectors at the center of the mainlobes, which by definition generatenon-clutter signals:(i) the local phase difference between the main and auxiliary receivebeams would be close to 0;(ii) the local magnitude ratio between the main and auxiliary receivebeams, wherein the auxiliary receive beam's signal is normalized so thatits maximal gain matches the maximal gain of the main receive beam,would be close to 1.0; and(iii) the local magnitude difference between the main and auxiliaryreceive beams, wherein the auxiliary receive beam's signal is normalizedso that its maximal gain matches the maximal gain of the main receivebeam, would be close to 0.0.

These assumptions may be derived from beamforming equations such as eq.(1) and eq. (2). A schematic illustration of two exemplary beampatterns, 51 and 52, defined in terms of gain as a function of theazimuth angle, is shown in FIG. 2. The resulting magnitude ratio betweenthe main and auxiliary receive beams, defined for a point-like reflectoras a function of the azimuth angle, can be seen in graph 55 in FIG. 3A.The resulting phase difference between the main and auxiliary receivebeams, defined for a point-like reflector as a function of the azimuthangle, can be seen in graph 56 in FIG. 3B. Note that the azimuth anglescale in FIG. 3B matches that of FIG. 3A, but does not match that ofFIG. 2.

Similarly, when using a main receive beam and an auxiliary receive,whose mainlobes do not point at the same spatial angle, one would expectthat for reflectors at the center of the mainlobes, which by definitiongenerate non-clutter signals:

(i) the local phase difference between the main and auxiliary receivebeams would be close to a constant, which may be calculated based on thespatial angles of the mainlobes of the main and auxiliary receive beams;and(ii) the local magnitude ratio between the main and auxiliary receivebeams, wherein the auxiliary receive beam's signal is normalized so thatits maximal gain matches the maximal gain of the main receive beam,would be close to a constant, which may be calculated based on thespatial angles of the mainlobes of the main and auxiliary receive beams.

It should be noted that any local phase difference and/or magnitudedifference and/or magnitude ratio may correspond to more than oneangular direction of the signal source. For example, when using a mainreceive beam and an auxiliary receive beam, whose mainlobes point at thesame spatial angle, wherein both the main receive beam and the auxiliaryreceive beam have a beam pattern which is symmetric about the center oftheir mainlobes (e.g., a conical beam pattern), at any given range withrespect to the phase center of the transducer array, one would expectthe local phase difference and magnitude difference and magnitude ratioto be equal along concentric circular curves, whose center falls at thecenterline of the mainlobes, wherein the curves lie in the correspondingrange with respect to the phase center of the transducer array. Somedeviations from this rule may be caused by physical effects within themedium, such as diffraction, which affect the actual beam pattern withinthe scanned region.

By using more than one of the following parameters: the local phasedifference, the local magnitude difference, and the local magnituderatio, one may reduce the ambiguity in the information regarding theangular direction of the signal source. By utilizing information frommultiple auxiliary receive beams for a given main receive beam, saidambiguity may be further reduced.

(b) For non-clutter signals, which originate from angular directionsapproximately corresponding by definition to the center of the mainreceive beam's mainlobe, one would expect the signal level and/or itsmagnitude as a function of spatial location, to change in a similarfashion in the main and auxiliary receive beams, assuming that themainlobes of the main and auxiliary receive beams sufficiently overlap,and that the data for the auxiliary beams is normalized so that theauxiliary beam's maximal gain would match the maximal gain of thecorresponding main beam. In that context, the signal level as a functionof spatial location may be defined, for example, by the local signallevel, the local spatial derivatives and/or the local normalized spatialderivatives, wherein a local normalized spatial derivative is defined asthe local spatial derivative divided by the local signal.

It should be noted that the lateral resolution of various receive beamsmay be different due to differences in the corresponding beam patterns,and that the axial resolution of various receive beams may be differentdue to dissimilar waveforms on transmission and/or reception. Therefore,in some embodiments, before comparing the signal level and/or itsmagnitude as a function of spatial location between different receivebeams, spatial filtering may be applied to one or more datasets, each ofwhich corresponding to a specific type of beam pattern and/or waveform,in order to make sure that the lateral and/or axial resolution of thelocal signal in the compared beams would match. When using multiplecommon configuration scanned data arrays, each of which corresponding toreceive beams with specific lateral and/or axial resolutions, spatialfiltering may be applied to one or more of the common configurationscanned data arrays before comparing the signal level and/or itsmagnitude as a function of spatial location between different receivebeams. For example, when using narrow main receive beams and wideauxiliary receive beams, wherein for each main receive beam there is acorresponding auxiliary beam whose mainlobe points at the samedirection, one may apply a spatial low-pass filter to the commonconfiguration scanned data arrays acquired with the main receive beam.

(c) Different target regions may have different dynamic models. Forexample, in echocardiography, the ribcage, which often causes clutterreflections, moves at much slower rates than the heart muscle within theregion of interest. Various temporal filters may be employed in order toestimate the signal levels for tissues with predefined dynamicproperties. The clutter suppression processes described herein may beapplied to specific frequency bands. Additionally or alternatively, onemay apply the clutter suppression processes described herein separatelyto outputs of several temporal filters, and then compound the results inorder to obtain the final dataset.

The compounding function parameters may be derived and/or utilized inone or more of the following manners (“compounding functionimplementation manners”):

(a) Computed for each entry into the scanned data array and each frameseparately, and applied accordingly.(b) Computed for a predefined subset of the frames, and in eachframe—for each entry into the scanned data array or group of adjacententries, and applied to all frames, wherein the parameters applied in aspecific frame may be the ones computed for the closest frame index, oralternatively, the closest frame index which is equal to or lower thanthe current frame index. Spatial interpolation and/or interpolation intime may also be used in order to derive the parameters applied in agiven frame.(c) Computed for a specific frame but for each entry into the scanneddata array or group of adjacent entries, and applied to all frames.(d) Computed for all frames but for a specific set of entries into thescanned data array, and applied to all entries into the scanned dataarray, wherein in each frame, the parameter value for each voxel isderived by interpolation and/or extrapolation, e.g., using nearestneighbor interpolation.(e) Computed for a predefined subset of the frames, whose length mayalso equal 1, and for a specific set of entries into the scanned dataarray, wherein in each frame, the parameter values are based on thosecomputed for the closest frame index, on those computed for the closestframe index which is equal to or lower than the current frame index, oron interpolation and/or extrapolation in time, and wherein spatialinterpolation and/or extrapolation is used for each entry into thescanned data array. The interpolation and/or extrapolation may utilizeany method known in the art.

Methods (b)-(e) are designed to reduce the time and/or memory complexityof the computations.

Generalized Metrics

In certain embodiments, one or more of the compounding functionparameters may be defined as adhering to at least one of the followingoptions:

(a) Indicative of the probability for the corresponding voxel to besubstantially affected by clutter effects only (i.e., almost all itsreceived energy originates from clutter effects), as derived usingcertain physical and/or physiologic assumptions.(b) Indicative of the probability for the corresponding voxel to besubstantially unaffected by clutter effects, as derived using certainphysical and/or physiologic assumptions. Note that, for each specificvoxel, the sum of the probabilities in (a) and (b) is always 1.0.(c) Indicative of the percentage of the received energy within thecorresponding voxel that originates from clutter effects.(d) Set to a certain constant, e.g., 0.0, if the corresponding voxel isnot significantly affected by clutter effects, and otherwise to adifferent constant, e.g., 1.0.

Such compounding function parameters are referred to herein as“generalized metrics.” They may be computed using any method known inthe art. In some embodiments, e.g., in option (d), segmentationtechniques may applied to at least one of the following (commonlyreferred to as the “segmentation technique parameters”): (a) the localreceived signal in one or more of the receive beams; and (b) the localvalue of one or more of the compounding function parameters.

For instance, classic edge detection (see, for example, U.S. Pat. No.6,716,175, by Geiser and Wilson, issued on Apr. 6, 2004, and titled“Autonomous boundary detection system for echocardiographic images”) orradial search techniques (see, for example, U.S. Pat. No. 5,457,754, byHan et al, issued on Oct. 10, 1995, and title “Method for automaticcontour extraction of a cardiac image”) may be used for segmentationpurposes. Such techniques may be combined with knowledge-basedalgorithms, aimed at performance enhancement, which may either beintroduced during post-processing, or as a cost-function, incorporatedwith the initial boundary estimation. Another example for an applicablesegmentation method is solving a constrained optimization problem, basedon active contour models (see, for example, a paper by Mishra et al.,entitled “A GA based approach for boundary detection of left ventriclewith echocardiographic image sequences,” Image and Vision Computing,vol. 21, 2003, pages 967-976, which is incorporated herein byreference). Additionally or alternatively, one may apply a predefinedthreshold to one or more of the segmentation technique parameters, inorder to obtain an initial division between voxels that aresignificantly affected by clutter effects and voxels that are not. Insome embodiments, iterative region growing techniques may then beemployed, wherein the predefined threshold for one or more of thesegmentation technique parameters may be changed between consecutiveiterations.

In further embodiments, the generalized metrics may be defined as realnumbers whose values range from 0.0 to 1.0, wherein the value 0.0 isassigned to voxels that are substantially unaffected by clutter effects,and the value 1.0 is assigned to voxels in which substantially all ofthe measured signal emanates from clutter effects (“normalizedgeneralized metrics”). In options (a)-(b) immediately above, thenormalized generalized metrics can be referred to as the estimatedprobability for the corresponding voxel to be substantially affected byclutter effects only; in option (c) immediately above, the normalizedgeneralized metrics can be referred to as the estimated percentage ofthe received energy within the corresponding voxel that originates fromclutter effects. Such definitions simplify signal calibration,especially when the compounding function is linear. In that context, theterm “signal calibration” refers to making sure that the signal energyin voxels completely unaffected by clutter effects changes as little aspossible due to the clutter suppression process, whereas the signalenergy in voxels whose sole energy source is clutter would be suppressedas much as possible by the clutter suppression process.

In certain embodiments, the generalized metric for a set of voxels in aset of frames may be derived using the following four-step process:

(a) Select one or more compounding function parameters to be used(“compounding function parameter set”). Some examples for compoundingfunction parameter sets:(i) parameter #1: the local magnitude ratio between two differentreceive beams (a main and an auxiliary beam);parameter #2: the local phase difference between two different receivebeams;(ii) parameter #1: the local ratio of spatial magnitude derivativesalong a first spatial axis (e.g., the X axis) between two differentreceive beams;parameter #2: the local ratio of spatial magnitude derivatives along asecond spatial axis (e.g., the Y axis) between two different receivebeams;(iii) parameter #1: the local ratio of normalized spatial magnitudederivatives along a first spatial axis (e.g., the X axis) between twodifferent receive beams;parameter #2: the local ratio of normalized spatial magnitudederivatives along a second spatial axis (e.g., the Y axis) between twodifferent receive beams;wherein a normalized local signal derivative is defined as the localsignal derivative divided by the local signal;(iv) parameter #1: the ratio of local magnitude spatial derivativesbetween two different receive beams, wherein the spatial derivative iscomputed using a “Laplacian filter”;parameter #2: the local magnitude ratio between two different beams;(v) parameter #1: the ratio of local normalized magnitude spatialderivatives between two different receive beams, wherein the normalizedspatial derivative is computed by dividing the local output of a“Laplacian filter” (applied to the local signal magnitude) by the localsignal magnitude;parameter #2: the local magnitude ratio between two different receivebeams; and(vi) parameter #1: the ratio of the normalized magnitude spans betweentwo different receive beams, wherein a normalized magnitude span isdefined as the difference between the maximal and minimal magnitudes(with or without outlier rejection) within a local region, divided bythe mean magnitude in said region, wherein the local region is definedso that its width is greater than its height;parameter #2: the ratio of the normalized magnitude spans between twodifferent beams, wherein the local region is defined so that its heightis greater than its width.(b) Define a model for at least one of:(i) the probability (or a simple function thereof, such as theprobability density function) for a voxel to be substantially affectedby clutter effects only;(ii) the probability (or a simple function thereof, such as theprobability density function) for a voxel to be substantially unaffectedby clutter effects; and(iii) the percentage of the energy within the corresponding voxel thatoriginates from clutter effects (or a simple function thereof);as a function of the compounding function parameter set values or thevalues for a subset of the compounding function parameter set(“compounding model”). In some cases, more than one compounding modelmay be employed, using one or more subsets of the compounding functionparameter set (“compounding function parameter subsets”). In thatcontext, the term “subset” also includes the option of the entire set.(c) Compute the compounding function parameter set for each relevantvoxel and/or its immediate vicinity. In that context, the immediatevicinity may be defined in space and/or in time over the scanned dataarray.(d) For each relevant voxel, use the one or more compounding models toestimate at least one of the following (or a simple function thereof):(i) the probability for its being substantially affected by cluttereffects only;(ii) the probability for its being substantially unaffected by cluttereffects; and(iii) the percentage of the energy within the corresponding voxel thatoriginates from clutter effects.

For each compounding model, the results for each relevant voxel may beused as a generalized metric or a normalized generalized metric for thatvoxel. Multiple generalized metrics may thus be defined for each voxel.

The compounding model may be defined using any method known in the art.Additionally or alternatively, it may be defined in one of more of thefollowing ways:

(a) A predefined model (“predefined compounding model”), which may betheoretically determined and/or experimentally derived. Such a model maybe based on the beam pattern of applicable transmit and/or receivebeams. For example, the beam width in one or more axes of each such beammay be taken into account. It may be independent of spatial location andtime, but it may also be dependent on spatial location and/or time. Forexample, one may use the fact that in the far-field, the nominaltransmit beam width changes as a function of the distance from theprobe, reaching a minimum at the focal length. Another example would beto use different models for different beam patterns, if more than onebeam pattern type is in use.(b) A simplified adaptive model (“simplified adaptive compoundingmodel”), based on an assumption regarding the shape of the probabilitydensity function (PDF) for a voxel to be substantially affected (orsubstantially unaffected) by clutter effects as a function of theapplicable compounding function parameter subset (“compounding modelPDF”). For example, a one-dimensional or multi-dimensional Gaussianfunction may be used, wherein the function's dimensions shouldcorrespond to the number of parameters in the compounding functionparameter subset.

The free parameters (in the example of the Gaussian function—the vectorof expectation values and the covariance matrix) may be determined bystatistical analysis of the values of the compounding function parameterset over some or all of the relevant voxels in the appropriate receivebeams. In some embodiments, all relevant voxels in all frames may beused at once. In other embodiments, the computations may be performedseparately for each frame and/or each group of voxels. In that context,groups of voxels may be separated in one or more of the scanned dataarray axes (e.g., different range swaths and/or different angularsections with respect to the broadside of the probe). In furtherembodiments, the computations may be performed for a subset of theframes and/or a subset of the voxels, wherein spatial and/or temporalinterpolation and/or extrapolation may be used to estimate thecompounding model PDF for any voxel. The interpolation and/orextrapolation may employ any method known in the art.

(c) An adaptive model (“adaptive compounding model”), wherein thecompounding model PDF, i.e., the PDF for a voxel to be substantiallyaffected (or substantially unaffected) by clutter effects as a functionof the applicable compounding function parameter subset, is computeddirectly, without any prior assumptions regarding the shape of the PDF.The compounding model PDF may be computed using any method known in theart. For example, a one- or multi-dimensional histogram of the values ofthe applicable compounding function parameter subset for all relevantvoxels may be calculated, wherein the histogram's dimensions shouldcorrespond to the number of parameters in the compounding functionparameter subset. In order to adhere to the definition of a PDF, thehistogram should be normalized such that the sum of all its elementswould equal 1.0. If one assumes that:(i) in most voxels, the primary source of the measured signal lies in aspatial angle close to the center of the mainlobe of the two-way beampattern (this is because the beam pattern provides higher gain to thisangular region, and most soft tissues have approximately the samereflective properties), and therefore voxels which are almost clutterfree are also more ubiquitous than voxels which are strongly affected byclutter; and(ii) the compounding function parameter subset provides informationregarding the local clutter levels and/or the angular direction of theprimary source of the measured signal;

the result may be used as an estimate of the PDF for a voxel to besubstantially unaffected by clutter effects. If one goes over eachhistogram element, and replaces its value by 1.0 minus the originalvalue, the result may be used as an estimate of the PDF for a voxel tobe substantially affected by clutter effects. Before computing thehistogram, one may apply outlier rejection in one or more of its axes,i.e., for one or more of the parameters in the compounding functionparameter subset, thus limiting the histogram's dynamic range. In someembodiments, the bins used in each histogram axis may be of a uniformwidth. In other embodiments, the bins used in one or more of thehistogram axes may be of adaptive width, wherein narrower bin widths areallotted to regions with higher PDF values. In even further embodiments,Voronoi diagrams may be used instead of histograms.

In certain embodiments, the compounding model PDF may be computed usingall relevant voxels in all frames at once. In other embodiments, thecompounding model PDF may be computed separately for each frame and/oreach group of voxels. In that context, groups of voxels may be separatedin one or more of the scanned data array axes (e.g., different rangeswaths and/or different angular sections with respect to the broadsideof the probe). In further embodiments, the compounding model PDF may becomputed for a subset of the frames and/or a subset of the voxels,wherein spatial and/or temporal interpolation and/or extrapolation maybe used to estimate the compounding model PDF for any voxel. Theinterpolation and/or extrapolation may utilize any method known in theart. For example, one may compute the compounding model PDF for thefirst frame out of every set of N_(f) consecutive frames, and use thatcompounding model PDF for all N_(f) consecutive frames.

Out of the three models, the predefined compounding model is expected tobe most computationally efficient, but it cannot take into accountmedium-specific artifacts, e.g., local attenuation, scattering,refraction and diffraction effects, which may occur along the beam,changing its effective beam pattern. The adaptive compounding model ismost accurate, yet it may be more computationally intensive. In order toreduce time complexity, one may use the compounding functionimplementation manners defined above.

In further embodiments, slight variations to the above describedcompounding models may allow one to produce normalized generalizedmetric values. As defined hereinabove, normalized generalized metricsare generalized metrics defined as real numbers whose value ranges from0.0 to 1.0, wherein the value 0.0 is assigned to voxels that aresubstantially unaffected by clutter effects, and the value 1.0 isassigned to voxels in which substantially all of the measured signalemanates from clutter effects. As a result, for predefined compoundingmodels, by correctly calibrating the predefined model one can obtainnormalized generalized metrics instead of un-normalized generalizedmetrics. For simplified adaptive compounding models or adaptivecompounding model, normalized generalized metrics may be obtainedinstead of un-normalized generalized metrics by utilizing theprobability for a voxel to be substantially affected by clutter effectsonly (as a function of the applicable compounding function parametersubset) instead of the compounding model PDF.

The probability for a voxel to be substantially affected by cluttereffects only (“compounding model probability function”) may be derivedfrom the compounding model PDF as follows:

Let {p_(n)} be the compounding function parameter subset, f_(PDF)({p_(n)}) be the compounding model PDF, and P({p_(n)}) be thecompounding model probability function.

For a compounding model PDF relating to the PDF for a voxel to besubstantially affected by clutter effects only, P({p_(n)}) for a given{p_(n)} may be estimated by the sum over all values of f_(PDF)({p_(n)})in the model, which are lower than or equal to the value of f_(PDF) forthe specific {p_(n)}. In cases where f_(PDF)({p_(n)}) is non-discrete,an integral should be used instead of a sum. Note that a PDF is bydefinition normalized so that its sum (or integral) over the entiremodel is 1.0, so that the compounding model probability function rangesfrom 0.0 to 1.0, as expected.

For a compounding model PDF relating to the PDF for a voxel to besubstantially unaffected by clutter effects, P({p_(n)}) for a given{p_(n)} may be estimated by the sum over all values of f_(PDF)({p_(n)})in the model, which are higher than or equal to the value of f_(PDF) forthe specific {p_(n)}. In cases where f_(PDF)({p_(n)}) is non-discrete,an integral should be used instead of a SUM.

Beam Pattern Measurements

In some embodiments, the compounding model PDF and/or the compoundingmodel probability function may be utilized for estimating the local beampattern within a given medium, and/or features of that beam pattern,such as the width of the mainlobe in one or more axes. For instance, ifthe transmit beam pattern is common to a main receive beam and itsassociated auxiliary receive beams, the receive beam pattern and/orfeatures of the receive beam pattern may be estimated. Appropriatemodels should be defined for that purpose, which may be theoreticallydetermined and/or experimentally derived (“beam pattern measurementmodels”).

For example, when using a single transmit beam, a single main receivebeam and a single auxiliary receive beam, wherein the azimuth beam widthof all beams is the same, but the elevation beam width of the mainreceive beam is lower than that of the transmit beam and the auxiliaryreceive beam, the compounding model PDF and/or the compounding modelprobability function are expected to contain information regarding theelevation beam width on receive. Similarly, when using a single transmitbeam, a single main receive beam and two auxiliary receive beams—onewider than the main receive beam in azimuth only, and the other inelevation only, one may use each auxiliary receive beam separately inorder to deduce information regarding the beam width on receive alongthe applicable axis.

In further embodiments, the generalized metric and/or the normalizedgeneralized metric may be utilized to estimate the angular direction ofthe primary energy source of the signal measured on receive for one ormore voxels in one or more frames, using beam pattern measurementmodels. For example, when using a single transmit beam, a single mainreceive beam and a single auxiliary receive beam, wherein the azimuthbeam width of all beams is the same, but the elevation beam width of themain receive beam is lower than that of the transmit beam and theauxiliary receive beam, the generalized metric and/or the normalizedgeneralized metric are expected to be indicative of the elevation anglefrom which most of the received energy originates for each voxel or eachframe. Various configurations of transmit beams, main receive beamsand/or auxiliary receive beams may be utilized for estimating thespatial angle of the main energy source for each voxel on receive.

Linear Compounding Functions

In some embodiments, the compounding function may be defined as a linearfunction of the local information for a main receive beam and one ormore associated auxiliary receive beams. For example, if a singleauxiliary receive beam is utilized, and a single normalized generalizedmetric m is calculated for each entry into the scanned data array f)and/or each frame n separately (or for a subset of the entries into thescanned data array and/or for a subset of the frames), the compoundingfunction may be defined as follows:

S _(out)({right arrow over (p)},n)=[1−m({right arrow over (p)},n)]S_(main)({right arrow over (p)},n)  (3)

wherein S_(main) is the local measured signal for the main receive beam,and S_(out) is the clutter suppressed local signal. This equationattenuates the signal for clutter affected voxels, nullifying the signalfor voxels in which substantially all of the measured signal emanatesfrom clutter effects, whereas voxels determined to be clutter free areunaffected.

In other embodiments, if a single auxiliary receive beam is utilized,and a single normalized generalized metric is calculated for each entryinto the scanned data array and each frame (or a subset thereof), thecompounding function may be a linear combination of the local signal forthe main receive beam, S_(main), and the auxiliary receive beam,S_(aux):

S _(out)({right arrow over (p)},n)=S _(main)({right arrow over(p)},n)−{tilde over (W)}({right arrow over (p)},n)S _(aux)({right arrowover (p)},n)  (4)

wherein {tilde over (W)} is a local complex weight (“clutter suppressionweight”), which may be calculated for each entry into the scanned dataarray and/or each frame (or a subset thereof). The clutter suppressionweight is computed so as to minimize the clutter contribution toS_(out), whereas the relevant signal contribution is retained as much aspossible.

If one normalizes the signal for the different auxiliary receive beamsS_(aux) _(i) , where i is the auxiliary beam index, so that the gain ofall auxiliary receive beams at the angular direction corresponding tothat of the center of the main receive beam would match that of the mainreceive beam, to obtain normalized auxiliary receive beams, {tilde over(S)}_(aux) _(i) , one can relate to the signal emanating from the centerof the main receive beam as the “clutter free” signal, S_(free), and theremaining signal for each receive beam may be referred to as the cluttersignal, denoted by S_(main) ^(clutter) for the main receive beam and{tilde over (S)}_(aux) _(i) ^(clutter) for the i'th normalized auxiliaryreceive beam:

S _(main)({right arrow over (p)},n)=S _(free)({right arrow over(p)},n)+S _(main) ^(clutter)({right arrow over (p)},n)

{tilde over (S)} _(aux) _(i) ({right arrow over (p)},n)=S _(free)({rightarrow over (p)},n)+{tilde over (S)} _(aux) _(i) ^(clutter)({right arrowover (p)},n)  (5)

The clutter suppression weight should therefore be computed in such away that S_(out) would be as close as possible to S_(free). Unlike themethod in eq. (3), which significantly attenuates the signal in voxelswhere most of the measured signal emanates from clutter, this method maybe able to extract relevant information from voxels in which clutter isthe main source of the measured signal.

In order to simplify the description below, let us rewrite eq. (4) usingnormalized auxiliary receive beam signals:

S _(out)({right arrow over (p)},n)=S _(main)({right arrow over(p)},n)−W({right arrow over (p)},n){tilde over (S)} _(aux)({right arrowover (p)},n)  (6)

wherein W is the local complex weight calculated using the normalizedauxiliary receive beam signals.

When using eq. (6), the magnitudes and/or phases of S_(out) may beskewed locally and/or globally. This is because the clutter suppressionweight aims at nullifying clutter contributions, but as an adverseeffect it also affects clutter free signals. For example, in a clutterfree voxel, S_(out)=(1−W)S_(free). In order to prevent these effects, anormalization factor may be introduces into eq. (6):

$\begin{matrix}{{S_{out}\left( {\overset{\rightarrow}{p},n} \right)} = \frac{{S_{main}\left( {\overset{\rightarrow}{p},n} \right)} - {{W\left( {\overset{\rightarrow}{p},n} \right)}{{\overset{\sim}{S}}_{aux}\left( {\overset{\rightarrow}{p},n} \right)}}}{1 - {W\left( {\overset{\rightarrow}{p},n} \right)}}} & (7)\end{matrix}$

Special attention should be paid to voxels for which W approximatelyequals 1.0, so that eq. (7) may yield very high signal magnitudes due to“divide-by-zero” effects. For example, S_(out) for such voxels may thedefined as equaling S_(main), so that the signal in such voxels remainsunchanged. Alternatively, S_(out) for such voxels may the set to themean or median signal level of S_(out) in the immediate vicinity of thecurrent voxel, wherein the immediate vicinity may be defined in spaceand/or in time, and wherein voxels for which the value provided by eq.(7) has been replaced are excluded from the computation. Such techniquesmay also be used in other occurrences of “divide-by-zero” effects.

In embodiments, the clutter suppression weight for each applicable entryinto the scanned data array {right arrow over (p)} and/or eachapplicable frame n may be estimated as follows:

(a) Define a region surrounding ({right arrow over (p)}, n), wherein theregion may be spatial (i.e., including a set of near-by entries into thescanned data array) and/or temporal (i.e., including a set ofconsecutive frames, whose index is close to n) (“spatial-temporalregion”).(b) Set the clutter suppression weight for ({right arrow over (p)}, n)to be a function of at least one of:(i) the values of one or more generalized metrics and/or normalizedgeneralized metrics within the spatial-temporal region; and(ii) the local complex ratio between the signal for the main receivebeam, S_(main), and the signal in one of the normalized auxiliaryreceive beams, {tilde over (S)}_(aux) _(i) within the spatial-temporalregion (“local complex main-auxiliary ratio”).

For example, the clutter suppression weight for ({right arrow over (p)},n) may be set to the local complex main-auxiliary ratio for to the voxelwithin the spatial-temporal region, whose generalized metric and/ornormalized generalized metric values are most indicative of significantclutter effects. This exemplary function is based on the followingassumptions:

(i) the clutter signal contribution changes slowly enough in spaceand/or in time, so that one may consider it to be approximately constantwithin small spatial-temporal regions; and(ii) speckle effects in ultrasonic imaging often cause local variationsin the clutter free signal levels even for approximately homogeneoustissue regions, so that one may expect at least one of the voxels in thespatial-temporal region to include low free signal levels.The voxel within the spatial-temporal region which is found to have themost substantial clutter effect may be treated as almost pure clutter,and therefore the local complex main-auxiliary ratio for that voxel isexpected to significantly reduce the local clutter level. For furtherclarification, if we substitute W({right arrow over (p)}, n)≈S_(main)^(clutter)({right arrow over (p)}, n)/{tilde over (S)}_(aux) _(i)^(clutter)({right arrow over (p)}, n) into eq. (7), we obtainS_(out)({right arrow over (p)}, n)≈S_(free)({right arrow over (p)}, n).(c) Special care should be taken to cases where the local complexmain-auxiliary ratio is very high due to “divide-by-zero” effects. Forexample, in voxels where the computed clutter suppression weight is veryhigh, one may replace it with a predefined value, e.g. 1.0 or 0.0. Inthe latter case, no clutter suppression will be applied to these voxels.

In further embodiments, one may compute a linear combination of theoutputs of a compounding function and the main receive beam signal,wherein the weights are based on the normalized generalized metric. Apotential benefit of this technique is reducing the effects of thecompounding function on clutter free voxels. For example, if we base ourmethod on eq. (7), we obtain:

$\begin{matrix}{{S_{out}\left( {\overset{\rightarrow}{p},n} \right)} = {{\left\lbrack {1 - {m\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack {S_{main}\left( {\overset{\rightarrow}{p},n} \right)}} + {{m\left( {\overset{\rightarrow}{p},n} \right)}\frac{{S_{main}\left( {\overset{\rightarrow}{p},n} \right)} - {{W\left( {\overset{\rightarrow}{p},n} \right)}{{\overset{\sim}{S}}_{aux}\left( {\overset{\rightarrow}{p},n} \right)}}}{1 - {W\left( {\overset{\rightarrow}{p},n} \right)}}}}} & (8)\end{matrix}$

In even further embodiments, two or more normalized generalized metricsmay be defined and utilized by the compounding function. For example,when using two normalized generalized metrics, m1 and m2, we can use thefollowing equation, wherein for each applicable entry into the scanneddata array and/or for each applicable frame, the sum of all coefficientsis, by definition, 1.0:

$\begin{matrix}{{S_{out}\left( {\overset{\rightarrow}{p},n} \right)} = {{{\left\lbrack {1 - {m_{1}\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack \left\lbrack {1 - {m_{2}\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack}{S_{main}\left( {\overset{\rightarrow}{p},n} \right)}} + {{{m_{1}\left( {\overset{\rightarrow}{p},n} \right)}\left\lbrack {1 - {m_{2}\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack}{S_{1}\left( {\overset{\rightarrow}{p},n} \right)}} + {\quad{{\left\lbrack {1 - {m_{1}\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack {m_{2}\left( {\overset{\rightarrow}{p},n} \right)}{S_{2}\left( {\overset{\rightarrow}{p},n} \right)}} + {{m_{1}\left( {\overset{\rightarrow}{p},n} \right)}{m_{2}\left( {\overset{\rightarrow}{p},n} \right)}{S_{1,2}\left( {\overset{\rightarrow}{p},n} \right)}}}}}} & (9)\end{matrix}$

wherein S₁ and S₂ are the outputs of clutter suppression processes basedon m1 and m2 respectively, applied using any method herein above; andS_(1,2) is the output of a clutter suppression process based on both m1and m2. This may be done, for example, by applying one of the aboveprocesses based on m1 to both the main receive beam and the auxiliaryreceive beams (i.e., each time a different beam is treated as the mainreceive beam), and then apply one of the above processes based on m2,substituting the receive beam signals in the appropriate equation withthe corresponding outputs of the clutter suppression process based onm1. Alternatively, S_(1,2) may be computed as a simple local statisticalfunction of S₁ and S₂, e.g., local minimum (between S₁ and S₂), localaverage (between S₁ and S₂), and so forth.

In other embodiments, two or more normalized generalized metrics may bedefined and utilized in the compounding function, wherein at least onenormalized generalized metric is used directly and at least onenormalized generalized metric is used to compute a clutter suppressionweight. For example, when using two normalized generalized metrics, m1and m2, wherein a clutter suppression weight W2 is computed using m2,the following equation may be used:

$\begin{matrix}{{S_{out}\left( {\overset{\rightarrow}{p},n} \right)} = {{{\left\lbrack {1 - {m_{1}\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack \left\lbrack {1 - {m_{2}\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack}{S_{main}\left( {\overset{\rightarrow}{p},n} \right)}} + {{m_{2}\left( {\overset{\rightarrow}{p},n} \right)}\frac{{S_{main}\left( {\overset{\rightarrow}{p},n} \right)} - {{W_{2}\left( {\overset{\rightarrow}{p},n} \right)}{{\overset{\sim}{S}}_{aux}\left( {\overset{\rightarrow}{p},n} \right)}}}{1 - {W_{2}\left( {\overset{\rightarrow}{p},n} \right)}}}}} & (10)\end{matrix}$

Note that compounding functions defined as linear combinations of thelocal signal for the main receive beam and the auxiliary receive beammay be thought of in at least one of two ways:(a) As algebraic expressions, designed to suppress clutter effects.(b) As adaptive beam-forming schemes, designed to introduce nulls intothe effective beam-pattern of the overall process, in directionscorresponding to those of clutter sources.

It should be understood that the mathematical equations are provided asdetailed examples only. Furthermore, even though the detailed examplesrefer to linear equations, non-linear equations may also be utilized.

Filtering of Generalized Metrics and/or Clutter Suppression Weights

Before utilizing a generalized metric and/or a normalized generalizedmetric for computing a clutter suppression weight and/or in acompounding function, it may undergo spatial filtering and/or temporalfiltering. Similarly, before utilizing a clutter suppression weight in acompounding function, it may also undergo spatial filtering and/ortemporal filtering. The filtering may be applied using the coordinatesystem of the scanned data array, wherein a time axis may be formed bycombining the scanned data arrays collected in multiple frames, whereinconsecutive frames are placed one alongside the other.

Any filter known in the art may be applied for that purpose. Forexample, a linear low-pass filter, band-pass filter of high-pass filtermay be applied in one or more of the spatial axes and/or the time axis.Additionally or alternatively, non-linear filters such as a medianfilter may be applied in one or more of the spatial axes and/or the timeaxis. Such filters are expected to reduce local irregularities in thegeneralized metrics, normalized generalized metrics and/or cluttersuppression weights.

Other non-linear filters, applied in one or more of the spatial axesand/or the time axis, may be considered. For example, one may use thelocal minimum, local maximum, local mean plus local standard deviation,local mean minus local standard deviation, a certain percentile of thelocal signal, and so on. When applied to a normalized generalizedmetric, a local minimum is expected to reduce misdetections of voxels asclutter affected, since the presence of voxels for which the normalizedgeneralized metric value is low (i.e., voxels determined to include lowclutter levels, and therefore mostly relevant information) also reducesthe normalized generalized metric value in their immediate vicinity.Similarly, when applied to a normalized generalized metric, a localmaximum is expected to increase the detection probability of clutteraffected voxels, since the presence of voxels for which the normalizedgeneralized metric value is high (i.e., voxels determined to includemostly clutter) also increases the normalized generalized metric valuein their immediate vicinity.

Further, slightly more complex, non-linear filters may also beconsidered. For example, for each local region in space and/or time, theoutput of the non-linear filter may be produced as follows: if the localmedian or the local mean or the local weighted mean is higher than apredefined value, use a certain statistical operator, andotherwise—another statistical operator. Examples for such statisticaloperators are the local minimum, local maximum, local mean, local meanplus local standard deviation, local mean minus local standarddeviation, a certain percentile of the local signal and so forth. Suchfilters, when applied to a normalized generalized metric, may allow one,for example, to increase the detection probability of clutter affectedvoxels in clutter affected regions and, at the same time, reduce theprobability of misdetections in clutter free regions, assuming thatclutter levels change relatively gradually in space and time.

Additionally or alternatively, a transfer function may be applied to ageneralized metric and/or a normalized generalized metric and/or aclutter suppression weight in various stages of the clutter suppressionprocess. Any transfer function known in the art may be used, withvarious purposes in mind. Transfer functions applied to a normalizedgeneralized metric, which by definition ranges from 0.0 to 1.0, shouldpreferably: (a) yield values between 0.0 and 1.0; and (b) for each pairof values, wherein the second value is higher than the first one, theresult of the transfer function for the second value should be higherthan or equal to the result of the transfer function for the firstvalue.

For example, a transfer function may be applied to the normalizedgeneralized metric m in eq. (3) and/or eq. (8) before using it forclutter suppression. Some exemplary transfer functions:

(a) Linear function with saturation; i.e., the transfer function islinear up to a predefined input value, and for higher input values—theoutput value is equal to that for the predefined input value.(b) Logarithmic function.(c) Exponential function.(d) Sigmoid-like function.

Post-Processing

Before utilizing the output of a compounding function for any furthercomputations, one may apply spatial and/or temporal processing to thatoutput. The processing may be applied using the coordinate system of thescanned data array, wherein a time axis may be formed by combining thescanned data arrays collected in multiple frames, wherein consecutiveframes are placed one alongside the other.

In some embodiments, the spatial and/or temporal processing may compriseof various spatial and/or temporal filters known in the art.

Additionally or alternatively, for each applicable entry into thescanned data array and/or each applicable frame, one may replace theoutput of the compounding function with the corresponding voxel of themain receive beam if a certain criterion is met. Alternatively, one mayreplace the local output of the compounding function with the meansignal, weighted mean signal or median signal within a small localregion of the compounding function output, defined in space and/or intime.

For example, the local value of the output of a compounding function maybe replaced by that of the corresponding voxel in the main receive beamif the local signal magnitude of the compounding function output isgreater than the local signal magnitude of the main receive beam. Such aprocess may be considered, because the clutter suppression process issupposed to reduce clutter levels, without amplifying clutter freeregions. For instance, high signal magnitudes may result from numericaleffects such as “divide-by-zero.”

Another exemplary post-process:

(a) Define a block within the scanned data array, centered orapproximately centered at the current voxel (“data block”). Depending onthe dimensions of the scanned data array (which may include spatial andpotentially also temporal axes), said block may be 1D, 2D, 3D orfour-dimensional (4D).(b) For each of the main receive beam and the output of the compoundingfunction, calculate the mean magnitude or the weighted mean magnitudewithin the data block, including or excluding the current voxel, anddivide the result by the magnitude of the current voxel (alternatively,the mean magnitude in a smaller block within the scanned data array maybe used).(c) The value of the current voxel within the compounding functionoutput should be replaced by that of the corresponding voxel in the mainreceive beam if the ratio between the output of the previous step forthe main receive beam and for the output of the compounding function ishigher-than (or lower-than) than a predefined constant. Alternatively,if the aforementioned criterion is met, the value of the current voxelwithin the compounding function output may be replaced by the meansignal, weighted mean signal or median signal within the data block ofthe compounding function output (either including or excluding thecurrent voxel).

Such a post process may be used, for instance, in order to take care ofsmall dark regions within the compounding function output, which appearin bright continuous (and larger) regions of the image.

Utilizing a Plurality of Main Receive Beams and/or a Plurality ofAuxiliary Receive Beams

In some embodiments, a plurality of main receive beams may be usedconcurrently, e.g., in MLA configurations. In such cases, one or moreauxiliary receive beams may be defined, wherein each auxiliary receivebeam is used in conjunction with one or more main receive beam.

In such cases, one possible clutter suppression method is to compute atleast one of the following “compounding process functions”:

(a) The compounding model PDF.(b) The generalized metrics and/or normalized generalized metrics.(c) The clutter suppression weights.(d) The compounding function.

for each applicable entry into the scanned data array and/or eachapplicable frame employing each pair of main receive beam and associatedauxiliary receive beam used (“main-auxiliary beam pair”) separately.

Additionally or alternatively, one may define compounding model PDFs (orother compounding process functions) which are common to more than onemain-auxiliary beam pair, e.g., for multiple main receive beams and onecommon auxiliary receive beam. In such cases, one may applytransformations to the data for each main-auxiliary beam pair within thegroup of pairs, which correspond to the estimated differences betweenvarious main-auxiliary beam pairs (“main-auxiliary beam pair datatransformations”). For example:

(a) When using a simplified adaptive compounding model and/or anadaptive compounding model, wherein a single compounding model PDF(“common compounding model PDF”) is used for multiple main-auxiliarybeam pairs, one may apply a tailored transformation to the data fromeach main-auxiliary beam pairs before using it to produce thecompounding model PDF (“main-auxiliary beam pair data transformation”).(b) When calculating the generalized metrics and/or normalizedgeneralized metrics, the common compounding model PDF may be transformedfor each main-auxiliary beam pair before using it for generalized metricand/or normalized generalized metric computation.

Examples for main-auxiliary beam pair data transformations, one or moreof which may be used when adding information to the common compoundingmodel PDF (similar but inverse transformation may be employed when usingthe common compounding model PDF values for generalized metric and/ornormalized generalized metric calculations):

(a) Before adding information based on a given main-auxiliary beam pairto the common compounding model PDF, a constant may be added to one ormore of the information's components, corresponding to one or more ofthe compounding function parameters referred to by the compounding modelPDF. This transformation is referred to as translation of thecompounding model PDF. For example, if two or more main receive beamsand one auxiliary receive beam are used, wherein the beam patterns ofthe main receive beams are identical except for spatial rotation aboutthe phase center of the transducer array, one may expect that thecompounding model PDF calculated for each of the two or moremain-auxiliary beam pairs would be similar, but translated in one ormore of the compounding function parameters referred to by thecompounding model PDF.(b) Before adding information based on a given main-auxiliary beam pairto the common compounding model PDF, either a linear or a non-lineartransformation may be applied to one or more of the information'scomponents, corresponding to one or more of the compounding functionparameters referred to by the compounding model PDF. The linear ornon-linear transformation may be applied to one or more compoundingfunction parameters separately, but it may also be applied to two ormore of the compounding function parameters combined. For example, iftwo or more main receive beams and one auxiliary receive beam are used,wherein the beam patterns of the mainlobes of all main receive beamspoint at the same direction, but their beam widths are different, onemay expect that the compounding model PDF calculated for each of the twoor more main-auxiliary beam pairs would be similar, butstretched/compressed along one or more of the compounding functionparameters referred to by the compounding model PDF.

In further embodiments, a plurality of auxiliary receive beams may beused concurrently, either with a single main receive beam or with aplurality of main receive beams. One or more of the following techniquesmay be employed is such cases:

(a) One or more of the aforementioned clutter suppression schemes may beapplied for each main-auxiliary beam pair separately. After that, foreach applicable small spatial volume and/or for each applicable frame, afunction may be applied to clutter suppression outputs for one or moreof the associated main-auxiliary beam pairs (“multi-pair combiningfunction”). Any function known in the art may be used, e.g., mean,maximum, minimum, median, a predefined percentile, weighted sum and soforth. For ordinal operators such as the median, the operatorcomputation should be based on sorting the relevant values according totheir magnitudes. The multi-pair combining function may also take intoaccount one or more generalized metrics and/or normalized generalizedmetrics, as computed for one or more of the main-auxiliary beam pairs,e.g., a weighted sum may be calculated, wherein higher weights areassigned to values for main-auxiliary beam pairs in which the signal isconsidered to include lower clutter levels.(b) For each applicable small spatial volume and/or for each applicableframe, one or more of the aforementioned clutter suppression schemes maybe applied for a specific main-auxiliary beam pair. The cluttersuppression output may then be referred to as a new (synthetic) main orauxiliary receive beam, to be used in conjunction with another auxiliaryand/or main receive beam as inputs to one or more of the above describedclutter suppression schemes. This process may be iteratively repeatedseveral times.(c) For each applicable small spatial volume and/or for each applicableframe, a compounding function may be applied to a main receive beams anda plurality of auxiliary receive beams, to produce a clutter suppressedoutput, wherein the compounding function employs one or more cluttersuppression weights. For example, if we extend eq. (6) to relate tomultiple auxiliary receive beams, we obtain:

S _(out)({right arrow over (p)},n)=S _(main)({right arrow over(p)},n)−Σ_(m) W _(m)({right arrow over (p)},n){tilde over (S)} _(aux)_(m) ({right arrow over (p)},n)  (11)

Similar equations may be written for any of the other cluttersuppression methods described herein.

The weight vector W_(m) in eq. (11), or similar coefficients in otherclutter suppression equations referring to multiple auxiliary receivebeam, may be estimated using any technique known in the art. Forexample, one may employ mean-square-error optimization, that is, for eq.(11), the weight vector W_(m) for each applicable small spatial and/ortemporal region may be set so as to minimize the expectation value ofthe output signal power, that is, E[|S_(main)({right arrow over (p)},n)−mWmp,nSauxmp,n2, wherein E is the expectation operator.

Another example would be to compute the weight vector W_(m) in eq. (11),or similar coefficients in other clutter suppression equations referringto multiple auxiliary receive beam, for each main-auxiliary beam pairseparately, without taking the other main-auxiliary beam pairs intoaccount. This methodology may be useful, for example, when the mainlobeof the main receive beam and the one or more auxiliary receive beamspoint at the same spatial direction, and for each main-auxiliary beampair the beam width of the main and auxiliary receive beams are the samealong a certain direction but different along another direction, whereinthe directions along with the beams have different widths in differentmain-auxiliary beam pairs may or may not be perpendicular. In suchcases, different main-auxiliary beam pairs may be utilized to suppressclutter originating from different spatial angles with respect to themainlobe of the main receive beam.

(d) For each applicable small spatial volume and/or for each applicableframe, a compounding function may be applied to a main receive beams anda plurality of auxiliary receive beams, to produce a clutter suppressedoutput, wherein the compounding function is linear, but does not employclutter suppression weights, for instance, as described in eq. (3). Inthis case, the generalized metric and/or normalized generalized metricmay be based on compounding function parameters which relate to morethan one auxiliary receive beam associated with one or more main receivebeams. For example, when using two auxiliary receive beams associatedwith a single main receive beam, one may use the following compoundingfunction parameters:(i) the local relative magnitude between the main receive beam and thefirst auxiliary receive beam;(ii) the local phase difference between the main receive beam and thefirst auxiliary receive beam;(iii) the local relative magnitude between the main receive beam and thesecond auxiliary receive beam; and(iv) the local phase difference between the main receive beam and thesecond auxiliary receive beam.(e) For each applicable small spatial volume and/or for each applicableframe, a compounding function may be applied to a main receive beams anda plurality of auxiliary receive beams, to produce a clutter suppressedoutput, wherein:(i) the compounding function employs one or more clutter suppressionweights; and(ii) one or more of the generalized metrics and/or normalizedgeneralized metrics may be based on compounding function parameterswhich relate to more than one auxiliary receive beam associated with oneor more main receive beams.(f) For each applicable small spatial volume and/or for each applicableframe, a compounding function may be applied to a main receive beams anda plurality of auxiliary receive beams, to produce a clutter suppressedoutput, wherein any compounding function may be employed.

As noted herein above, one or more receive beams may be employed as bothmain receive beams and as auxiliary receive beams, within the context ofdifferent main-auxiliary beam pairs.

Utilizing Main and Auxiliary Receive Beams Having Different PhaseCenters

In some embodiments, for each main receive beam, the phase centers ofall associated auxiliary receive beams are identical to that of thecorresponding main receive beam. In other embodiments, there is at leastone main-auxiliary beam pair for which the main and auxiliary receivebeams have different phase centers (“non-concentric main-auxiliary beampair”).

An exemplary use for such configurations is suppression of reverberationeffects, in addition to sidelobe clutter effects. Reverberationartifacts measured for different receive phase centers are expected tobe different, so that compounding functions may suppress both sidelobeclutter and reverberations.

In non-concentric main-auxiliary beam pairs, the scanned data arrayassociated with the auxiliary receive beam may be translated and/orrotated with respect to the scanned data array associated with thecorresponding main receive beam. In certain embodiments, spatialregistration may be applied between the scanned data arrays associatedwith the main and auxiliary receive beams in non-concentricmain-auxiliary beam pair (“main-auxiliary beam pair registration”). Theregistration may be performed for the entire cine-loop, for each frameseparately, in each frame—for each entry into the scanned data array orgroup of adjacent entries, or for each entry into the scanned data arrayor group of adjacent entries in the entire cine-loop. The spatialregistration may be performed using any technique known in the art. Insome cases, the spatial registration may be rigid, accounting for globaltranslations and/or global rotation. Additionally or alternatively, thespatial registration may be non-rigid, also taking into account localdeformations resulting from physical effects within the medium. Incertain cases, some translation parameters, such as the globaltranslation and the global rotation, may be directly derived from therelative geometry of the centerlines of the main and auxiliary receivebeams, which are known and controllable. Additionally or alternatively,some translation parameters may be initially estimated according to therelative geometry of the centerlines of the main and auxiliary receivebeams, and then fine-tuned according to the respective acquired scanneddata arrays.

In non-concentric main-auxiliary beam pairs, the spatial angle betweenthe center of a voxel and the centerline of the main receive beam (“mainspatial angle”) is typically different than the spatial angle betweenthe center of said voxel and the centerline of the respective auxiliaryreceive beam (“auxiliary spatial angle”). Moreover, the differencebetween the main spatial angle and the auxiliary spatial angle may bedifferent for different voxels. As a result, the local relativemagnitude and/or the local phase difference between the main and theauxiliary receive beam, as calculated for reflectors along thecenterline of the main receive beam, may change as a function of range(“parameter range variability”). In certain embodiments, the compoundingfunction may be adjusted in order to handle this parameter rangevariability. For example, one may utilize main-auxiliary beam pair datatransformations, as defined hereinabove, wherein the transformation forsome voxels may be derived from the main spatial angle and auxiliaryspatial angle for each voxel or group of voxels, wherein the mainspatial angles and auxiliary spatial angles may be computed based on thegeometric location of the center of the voxel with respect to the phasecenters of the main and auxiliary receive beams, and/or from the outputsof the main-auxiliary beam pair registration process, if applied.

Clutter Suppression Based on Both Clutter De-Correlation Times andSpatial-Temporal Analysis

U.S. Pat. No. 8,045,777, also referenced herein above, describes aclutter suppression method based on computing a linear combination ofthe signals in a main receive beam and in one or more auxiliary receivebeams, wherein the coefficients of the linear combination (“temporalclutter suppression weights”) are computed so as to suppress reflectionsdue to clutter, wherein a reflection is determined as associated withclutter if the determined de-correlation time is above a specifiedthreshold (“temporal clutter suppression process”). When a singleauxiliary receive beam is associated with each main receive beam, thefollowing equation may be used:

S _(out)({right arrow over (p)},n)=S _(main)({right arrow over (p)},n)−W_(t)({right arrow over (p)},n)n){tilde over (S)} _(aux)({right arrowover (p)},n)  (12)

wherein the complex weight W_(t) for each applicable entry into thescanned data array {right arrow over (p)} and/or each applicable frame nis computed as follows:

$\begin{matrix}{{W_{t}\left( {\overset{\rightarrow}{p},n} \right)} = \frac{F_{t}\left\lbrack {S_{main}\left( {\overset{\rightarrow}{p},n} \right)} \right\rbrack}{F_{t}\left\lbrack {{\overset{\sim}{S}}_{aux}\left( {\overset{\rightarrow}{p},n} \right)} \right\rbrack}} & (13)\end{matrix}$

where Ft is a temporal low-pass filter, applied for each applicableentry into the scanned data array and/or each applicable frame.

This technique may be combined with the clutter suppression techniqueshereof, using compounding functions which also depend on temporalclutter suppression weights. For example, one or more generalizedmetrics and/or normalized generalized metrics may be used in addition toone or more temporal clutter suppression weights, as can be seen in thefollowing exemplary equation:

$\begin{matrix}{{S_{out}\left( {\overset{\rightarrow}{p},n} \right)} = {{\left\lbrack {1 - {m\left( {\overset{\rightarrow}{p},n} \right)}} \right\rbrack {S_{main}\left( {\overset{\rightarrow}{p},n} \right)}} + {{m\left( {\overset{\rightarrow}{p},n} \right)}\frac{{S_{main}\left( {\overset{\rightarrow}{p},n} \right)} - {{W_{t}\left( {\overset{\rightarrow}{p},n} \right)}{{\overset{\sim}{S}}_{aux}\left( {\overset{\rightarrow}{p},n} \right)}}}{1 - {W_{t}\left( {\overset{\rightarrow}{p},n} \right)}}}}} & (14)\end{matrix}$

wherein m is a normalized generalized metric, used so as to adjust theeffect of the temporal clutter suppression process in each voxel basedon the estimated local clutter contribution.

1. A method of ultrasound imaging, said method comprising: transmittingultrasound radiation towards a target and receiving reflections of theultrasound radiation from a region of the target in a main reflectedsignal and one or more auxiliary reflected signals, wherein each one ofthe reflected signals comprises an input dataset, the main reflectedsignal being associated with a main reflected beam and the one or moreauxiliary reflected signals are associated with one or more auxiliaryreflected beams, the main reflected beam and the auxiliary reflectedbeams being pointed in the same direction and having different anddistinct beam patterns; compounding the input datasets from the mainreflected signal and one or more auxiliary reflected signals, by the useof a compounding function, said compounding function using parametersderived from spatial analysis of the input datasets, said compoundingfunction parameters being derived from at least one of: (a) A localphase difference and/or magnitude difference and/or magnitude ratioand/or complex signal ratio between different receive beams, and/orfunctions of one or more of the aforementioned parameters; (b) Spatialfunctions of the local phase difference and/or magnitude differenceand/or magnitude ratio and/or complex signal ratio between differentbeams; (c) A local difference and/or ratio between spatial functions ofthe local magnitude and/or phase and/or complex signal in the differentreceive beams; (d) Applying multiple temporal band-pass filters to thelocal magnitude and/or phase and/or complex signal in the differentreceive beams, applying spatial analysis to the outputs of multipletemporal band-pass filters and compounding the results; (e) Applyingmultiple temporal band-pass filters to the local phase difference and/ormagnitude difference and/or magnitude ratio and/or complex signal ratiobetween different receive beams, applying spatial analysis to theoutputs of multiple temporal band-pass filters and compounding theresults; or (f) Spatial-temporal functions of the local phase differenceand/or magnitude difference and/or magnitude ratio between differentbeams, wherein spatial derivatives of the local phase difference and/ormagnitude difference and/or magnitude ratio between different beams, areapplied after local temporal filtering; or wherein the local differenceand/or ratio between spatial functions of the magnitude and/or phaseand/or complex signal in the different beams are applied after localtemporal filtering.
 2. The method of claim 1 wherein the methodcomprises receiving a plurality of main reflected signals and associatedauxiliary signals.
 3. The method of claim 1 wherein the main reflectedsignal and the auxiliary reflected signal have the same centerfrequency.
 4. The method of claim 1, wherein the output of thecompounding function for any region of the target in anytime frame maybe used instead of the measured main reflected signal in any furtherprocessing performed by ultrasound scanners, said further processingbeing at least one of: (a) Log-compression; (b) Time-gain control; (c)Global gain control; or (d) Doppler processing, wherein said Dopplerprocessing is one of: color-Doppler flow imaging; tissue-Dopplerimaging; pulse-Doppler studies; and/or continuous-wave Doppler studies.5. The method of claim 1, wherein the main reflected signal and the oneor more auxiliary reflected signals result from a single transmittedpulse or from two or more transmitted pulses.
 6. The method of claim 1,wherein the main reflected signal and the one or more auxiliaryreflected signals are obtained in one of the following ways: (a) Allobtained concurrently, wherein the main reflected signal and the one ormore auxiliary reflected signals are all associated with the sametransmit pulse; (b) All obtained at different times, wherein each of themain reflected signal and the one or more auxiliary reflected signals isassociated with a different transmit pulse; or (c) Divided into groupsof reflected signals, wherein each group of reflected signal maycomprise a main reflected signal and/or one or more auxiliary reflectedsignals, wherein each group of reflected signals is obtainedconcurrently but different groups are obtained at different times, andwherein each group of reflected signals is associated with a differenttransmit pulse.
 7. The method of claim 1, wherein the local phasedifference and/or magnitude difference and/or magnitude ratio areindicative of the angular direction of the signal source.
 8. The methodof claim 1, wherein for non-clutter signals which originate from angulardirections approximately corresponding to the center of the mainlobe ofthe beam associated with the main reflected signal, the signal leveland/or the magnitude of the signal level change as a function of spatiallocation approximately the same way in the main reflected signal and inthe one or more auxiliary reflected signals, wherein the change as afunction of spatial location is defined by at least one of: the localsignal level; the local spatial derivatives; and the local normalizedspatial derivatives, said local normalized spatial derivative beingdefined by the local spatial derivative divided by the local signal. 9.The method of claim 1, wherein temporal filters are employed to estimatethe signal levels for target elements with predefined dynamicproperties.
 10. The method of claim 1, wherein the dataset collected byone or more reflected signals, in each time-swath or frame, is organizedin a 1D, 2D or 3D scanned data array and the compounding functionparameter is computed for each entry into the scanned data array andeach frame, separately.
 11. The method of claim 1, wherein the datasetcollected by one or more reflected signals, in each time-swath or frame,is organized in a 1D, 2D or 3D scanned data array and the compoundingfunction parameter is computed for a predefined subset of the frames,and in each frame, for each entry into the scanned data array or groupof adjacent entries, and applied to all frames, wherein the parametersapplied in a specific frame are the ones computed for the closest frameindex, or alternatively, the closest frame index which is equal to orlower than the current frame index.
 12. The method of claim 11, whereinspatial interpolation and/or interpolation in time is used in order toderive the compounding function parameters applied in a given frame. 13.The method of claim 1, wherein the dataset collected by one or morereflected signals, in each time-swath or frame, is organized in a 1D, 2Dor 3D scanned data array and the compounding function parameter iscomputed for a specific frame and for each entry into the scanned dataarray or group of adjacent entries, and applied to all frames.
 14. Themethod of claim 1, wherein the dataset collected by one or morereflected signals, in each time-swath or frame, is organized in a 1D, 2Dor 3D scanned data array and the compounding function parameter iscomputed for all frames and for a specific set of entries into thescanned data array, and applied to all entries into the scanned dataarray, wherein in each frame, the parameter value for each voxel isderived by interpolation and/or extrapolation.
 15. The method of claim1, wherein the dataset collected by one or more reflected signals, ineach time-swath or frame, is organized in a 1D, 2D or 3D scanned dataarray and the compounding function parameter is computed for apredefined subset of the frames, and for a specific set of entries intothe scanned data array, wherein in each frame, the parameter values arebased on those computed for the closest frame index, on those computedfor the closest frame index which is equal to or lower than the currentframe index, or on interpolation and/or extrapolation in time, andwherein spatial interpolation and/or extrapolation is used for eachentry into the scanned data array.
 16. The method of claim 1, whereinthe compounding function parameter is a generalized metric and isindicative of the probability for the corresponding voxel to besubstantially affected by clutter effects only, as derived using certainphysical and/or physiologic assumptions.
 17. The method of claim 1,wherein the compounding function parameter is a generalized metric andis indicative of the probability for the corresponding voxel to besubstantially unaffected by clutter effects, as derived using certainphysical and/or physiologic assumptions.
 18. The method of claim 1,wherein the compounding function parameter is a generalized metric andis indicative of the percentage of the received energy within thecorresponding voxel that originates from clutter effects.
 19. The methodof claim 1, wherein the compounding function parameter is a generalizedmetric and is set to a certain constant, if the corresponding voxel isnot significantly affected by clutter effects, and otherwise to adifferent constant.
 20. The method of claim 19, wherein segmentationtechniques are applied to the local reflected signal in one or more ofthe main reflected signal and the one or more auxiliary reflectedsignals.
 21. The method of claim 19, wherein segmentation techniques areapplied to the local value of one or more of the compounding functionparameters.
 22. The method of claim 16, wherein the generalized metricis a normalized generalized metric and is defined as real numbers whosevalues range from 0.0 to 1.0, wherein the value 0.0 is assigned tovoxels that are substantially unaffected by clutter effects, and thevalue 1.0 is assigned to voxels in which substantially all of themeasured signal emanates from clutter effects, and wherein thenormalized generalized metric is the estimated probability for thecorresponding voxel to be substantially affected by clutter effectsonly.
 23. The method of claim 17, wherein the generalized metric is anormalized generalized metric and is defined as real numbers whosevalues range from 0.0 to 1.0, wherein the value 1.0 is assigned tovoxels that are substantially unaffected by clutter effects, and thevalue 0.0 is assigned to voxels in which substantially all of themeasured signal emanates from clutter effects, and wherein thenormalized generalized metric is the estimated probability for thecorresponding voxel to be substantially unaffected by clutter effects.24. The method of claim 18, wherein the generalized metric is anormalized generalized metric and is defined as real numbers whosevalues range from 0.0 to 1.0, wherein one of the following applies: (a)The value 0.0 is assigned to voxels that are substantially unaffected byclutter effects, and the value 1.0 is assigned to voxels in whichsubstantially all of the measured signal emanates from clutter effects,and wherein the normalized generalized metric is the estimatedpercentage of the received energy within the corresponding voxel thatoriginates from clutter effects, divided by 100; or (b) The value 1.0 isassigned to voxels that are substantially unaffected by clutter effects,and the value 0.0 is assigned to voxels in which substantially all ofthe measured signal emanates from clutter effects, and wherein thenormalized generalized metric is 100 minus the estimated percentage ofthe received energy within the corresponding voxel that originates fromclutter effects, divided by
 100. 25. The method of claim 19, wherein thegeneralized metric is a normalized generalized metric and is defined asreal numbers whose values range from 0.0 to 1.0, wherein one of thefollowing applies: (a) The value 0.0 is assigned to voxels that aresubstantially unaffected by clutter effects, and the value 1.0 isassigned to voxels in which substantially all of the measured signalemanates from clutter effects; or (b) The value 1.0 is assigned tovoxels that are substantially unaffected by clutter effects, and thevalue 0.0 is assigned to voxels in which substantially all of themeasured signal emanates from clutter effects.
 26. The method of claim16, wherein the generalized metric for a set of voxels in a set offrames is derived using the following four-step process: (a) Select oneor more compounding function parameters to be used (compounding functionparameter set) from the following parameters: (i) parameter #1: thelocal magnitude ratio between two different reflected beams (a main andan auxiliary beam); parameter #2: the local phase difference between twodifferent reflected beams; (ii) parameter #1: the local ratio of spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of spatial magnitudederivatives along a second spatial axis between two different reflectedbeams; (iii) parameter #1: the local ratio of normalized spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of normalized spatialmagnitude derivatives along a second spatial axis between two differentreflected beams; wherein a normalized local signal derivative is definedas the local signal derivative divided by the local signal; (iv)parameter #1: the ratio of local magnitude spatial derivatives betweentwo different receive beams, wherein the spatial derivative is computedusing a “Laplacian filter”; parameter #2: the local magnitude ratiobetween two different beams; (v) parameter #1: the ratio of localnormalized magnitude spatial derivatives between two different reflectedbeams, wherein the normalized spatial derivative is computed by dividingthe local output of a “Laplacian filter” (applied to the local signalmagnitude) by the local signal magnitude; parameter #2: the localmagnitude ratio between two different reflected beams; and (vi)parameter #1: the ratio of the normalized magnitude spans between twodifferent reflected beams, wherein a normalized magnitude span isdefined as the difference between the maximal and minimal magnitudes(with or without outlier rejection) within a local region, divided bythe mean magnitude in said region, wherein the local region is definedso that its width is greater than its height; parameter #2: the ratio ofthe normalized magnitude spans between two different beams, wherein thelocal region is defined so that its height is greater than its width.(b) Define a model for at least one of: (i) the probability (or a simplefunction thereof, such as the probability density function) for a voxelto be substantially affected by clutter effects only; (ii) theprobability (or a simple function thereof, such as the probabilitydensity function) for a voxel to be substantially unaffected by cluttereffects; and (iii) the percentage of the energy within the correspondingvoxel that originates from clutter effects (or a simple functionthereof); as a function of the compounding function parameter set valuesor the values for a subset of the compounding function parameter set(compounding model); (c) Compute the compounding function parameter setfor each relevant voxel and/or its immediate vicinity, said immediatevicinity being defined in space and/or in time; (d) For each relevantvoxel, use the one or more compounding models to estimate at least oneof the following (or a simple function thereof): (i) the probability forits being substantially affected by clutter effects only; (ii) theprobability for its being substantially unaffected by clutter effects;and (iii) the percentage of the energy within the corresponding voxelthat originates from clutter effects.
 27. The method of claim 17,wherein the generalized metric for a set of voxels in a set of frames isderived using the following four-step process: (a) Select one or morecompounding function parameters to be used (compounding functionparameter set) from the following parameters: (i) parameter #1: thelocal magnitude ratio between two different reflected beams (a main andan auxiliary beam); parameter #2: the local phase difference between twodifferent reflected beams; (ii) parameter #1: the local ratio of spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of spatial magnitudederivatives along a second spatial axis between two different reflectedbeams; (iii) parameter #1: the local ratio of normalized spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of normalized spatialmagnitude derivatives along a second spatial axis between two differentreflected beams; wherein a normalized local signal derivative is definedas the local signal derivative divided by the local signal; (iv)parameter #1: the ratio of local magnitude spatial derivatives betweentwo different receive beams, wherein the spatial derivative is computedusing a “Laplacian filter”; parameter #2: the local magnitude ratiobetween two different beams; (v) parameter #1: the ratio of localnormalized magnitude spatial derivatives between two different reflectedbeams, wherein the normalized spatial derivative is computed by dividingthe local output of a “Laplacian filter” (applied to the local signalmagnitude) by the local signal magnitude; parameter #2: the localmagnitude ratio between two different reflected beams; and (vi)parameter #1: the ratio of the normalized magnitude spans between twodifferent reflected beams, wherein a normalized magnitude span isdefined as the difference between the maximal and minimal magnitudes(with or without outlier rejection) within a local region, divided bythe mean magnitude in said region, wherein the local region is definedso that its width is greater than its height; parameter #2: the ratio ofthe normalized magnitude spans between two different beams, wherein thelocal region is defined so that its height is greater than its width.(b) Define a model for at least one of: (i) the probability (or a simplefunction thereof, such as the probability density function) for a voxelto be substantially affected by clutter effects only; (ii) theprobability (or a simple function thereof, such as the probabilitydensity function) for a voxel to be substantially unaffected by cluttereffects; and (iii) the percentage of the energy within the correspondingvoxel that originates from clutter effects (or a simple functionthereof); as a function of the compounding function parameter set valuesor the values for a subset of the compounding function parameter set(compounding model); (c) Compute the compounding function parameter setfor each relevant voxel and/or its immediate vicinity, said immediatevicinity being defined in space and/or in time; (d) For each relevantvoxel, use the one or more compounding models to estimate at least oneof the following (or a simple function thereof): (i) the probability forits being substantially affected by clutter effects only; (ii) theprobability for its being substantially unaffected by clutter effects;and (iii) the percentage of the energy within the corresponding voxelthat originates from clutter effects.
 28. The method of claim 18,wherein the generalized metric for a set of voxels in a set of frames isderived using the following four-step process: (a) Select one or morecompounding function parameters to be used (compounding functionparameter set) from the following parameters: (i) parameter #1: thelocal magnitude ratio between two different reflected beams (a main andan auxiliary beam); parameter #2: the local phase difference between twodifferent reflected beams; (ii) parameter #1: the local ratio of spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of spatial magnitudederivatives along a second spatial axis between two different reflectedbeams; (iii) parameter #1: the local ratio of normalized spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of normalized spatialmagnitude derivatives along a second spatial axis between two differentreflected beams; wherein a normalized local signal derivative is definedas the local signal derivative divided by the local signal; (iv)parameter #1: the ratio of local magnitude spatial derivatives betweentwo different receive beams, wherein the spatial derivative is computedusing a “Laplacian filter”; parameter #2: the local magnitude ratiobetween two different beams; (v) parameter #1: the ratio of localnormalized magnitude spatial derivatives between two different reflectedbeams, wherein the normalized spatial derivative is computed by dividingthe local output of a “Laplacian filter” (applied to the local signalmagnitude) by the local signal magnitude; parameter #2: the localmagnitude ratio between two different reflected beams; and (vi)parameter #1: the ratio of the normalized magnitude spans between twodifferent reflected beams, wherein a normalized magnitude span isdefined as the difference between the maximal and minimal magnitudes(with or without outlier rejection) within a local region, divided bythe mean magnitude in said region, wherein the local region is definedso that its width is greater than its height; parameter #2: the ratio ofthe normalized magnitude spans between two different beams, wherein thelocal region is defined so that its height is greater than its width.(b) Define a model for at least one of: (i) the probability (or a simplefunction thereof, such as the probability density function) for a voxelto be substantially affected by clutter effects only; (ii) theprobability (or a simple function thereof, such as the probabilitydensity function) for a voxel to be substantially unaffected by cluttereffects; and (iii) the percentage of the energy within the correspondingvoxel that originates from clutter effects (or a simple functionthereof); as a function of the compounding function parameter set valuesor the values for a subset of the compounding function parameter set(compounding model); (c) Compute the compounding function parameter setfor each relevant voxel and/or its immediate vicinity, said immediatevicinity being defined in space and/or in time; (d) For each relevantvoxel, use the one or more compounding models to estimate at least oneof the following (or a simple function thereof): (i) the probability forits being substantially affected by clutter effects only; (ii) theprobability for its being substantially unaffected by clutter effects;and (iii) the percentage of the energy within the corresponding voxelthat originates from clutter effects.
 29. The method of claim 19,wherein the generalized metric for a set of voxels in a set of frames isderived using the following four-step process: (a) Select one or morecompounding function parameters to be used (compounding functionparameter set) from the following parameters: (i) parameter #1: thelocal magnitude ratio between two different reflected beams (a main andan auxiliary beam); parameter #2: the local phase difference between twodifferent reflected beams; (ii) parameter #1: the local ratio of spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of spatial magnitudederivatives along a second spatial axis between two different reflectedbeams; (iii) parameter #1: the local ratio of normalized spatialmagnitude derivatives along a first spatial axis between two differentreflected beams; parameter #2: the local ratio of normalized spatialmagnitude derivatives along a second spatial axis between two differentreflected beams; wherein a normalized local signal derivative is definedas the local signal derivative divided by the local signal; (iv)parameter #1: the ratio of local magnitude spatial derivatives betweentwo different receive beams, wherein the spatial derivative is computedusing a “Laplacian filter”; parameter #2: the local magnitude ratiobetween two different beams; (v) parameter #1: the ratio of localnormalized magnitude spatial derivatives between two different reflectedbeams, wherein the normalized spatial derivative is computed by dividingthe local output of a “Laplacian filter” (applied to the local signalmagnitude) by the local signal magnitude; parameter #2: the localmagnitude ratio between two different reflected beams; and (vi)parameter #1: the ratio of the normalized magnitude spans between twodifferent reflected beams, wherein a normalized magnitude span isdefined as the difference between the maximal and minimal magnitudes(with or without outlier rejection) within a local region, divided bythe mean magnitude in said region, wherein the local region is definedso that its width is greater than its height; parameter #2: the ratio ofthe normalized magnitude spans between two different beams, wherein thelocal region is defined so that its height is greater than its width.(b) Define a model for at least one of: (i) the probability (or a simplefunction thereof, such as the probability density function) for a voxelto be substantially affected by clutter effects only; (ii) theprobability (or a simple function thereof, such as the probabilitydensity function) for a voxel to be substantially unaffected by cluttereffects; and (iii) the percentage of the energy within the correspondingvoxel that originates from clutter effects (or a simple functionthereof); as a function of the compounding function parameter set valuesor the values for a subset of the compounding function parameter set(compounding model); (c) Compute the compounding function parameter setfor each relevant voxel and/or its immediate vicinity, said immediatevicinity being defined in space and/or in time; (d) For each relevantvoxel, use the one or more compounding models to estimate at least oneof the following (or a simple function thereof): (i) the probability forits being substantially affected by clutter effects only; (ii) theprobability for its being substantially unaffected by clutter effects;and (iii) the percentage of the energy within the corresponding voxelthat originates from clutter effects.
 30. The method of claim 1, whereinthe compounding function parameter is a generalized metric, wherein acompounding function parameter set is a set of one or more compoundingfunction parameters to be used, wherein a compounding function parametersubset is either a subset or a full set of a compounding functionparameter set, wherein a compounding model is a model of at least oneof: (i) the probability (or a simple function thereof, such as theprobability density function) for a voxel to be substantially affectedby clutter effects only; (ii) the probability (or a simple functionthereof, such as the probability density function) for a voxel to besubstantially unaffected by clutter effects; and (iii) the percentage ofthe energy within the corresponding voxel that originates from cluttereffects (or a simple function thereof); as a function of the compoundingfunction parameter set values or the values for a subset of thecompounding function parameter set, and wherein the compounding modelmay be defined in one of more of the following ways: (a) A predefinedmodel, which may be theoretically determined and/or experimentallyderived; (b) A simplified adaptive model, based on an assumptionregarding the shape of the probability density function for a voxel tobe substantially affected (or substantially unaffected) by cluttereffects as a function of the applicable compounding function parametersubset (compounding model probability density function); (c) An adaptivemodel, wherein the compounding model probability density function iscomputed directly, without any prior assumptions regarding the shape ofthe probability density function.
 31. The method of claim 30, whereinvoxels that are substantially unaffected by clutter are assumed to bemore ubiquitous than voxels that are strongly affected by clutter. 32.The method of claim 30, wherein the compounding function parametersubset provides information regarding the local clutter levels and/orthe angular direction of the primary source of the measured signalenergy.
 33. The method of claim 30, wherein the compounding modelprobability density function is computed using all relevant voxels inall frames at once.
 34. The method of claim 30, wherein the datasetcollected by one or more reflected signals, in each time-swath or frame,is organized in a 1D, 2D or 3D scanned data array, and wherein thecompounding model probability density function is computed separatelyfor each frame and/or each group of voxels, wherein the groups of voxelsare separated in one or more scanned data array axes.
 35. The method ofclaim 30, wherein the compounding model probability density function iscomputed for a subset of the frames and/or a subset of the voxels,wherein spatial and/or temporal interpolation and/or extrapolation isused to estimate the compounding model probability density function forany voxel.
 36. The method of claim 30, wherein the compounding modelprobability density function is utilized for estimating the local beampattern within a given target region, and/or features of that beampattern.
 37. The method of claim 22, wherein the compounding function isdefined as a linear function of the local information for a mainreflected beam and one or more associated auxiliary beams and whereinthe compounding function is defined by the following equation:S _(out)=[1−m]S _(main) wherein m is the normalized generalized metric;S_(main) is the local measured signal for the main reflected beam; andS_(out) is the clutter suppressed local signal.
 38. The method of claim1, wherein the compounding function is defined by a linear combinationof the local signal for the main reflected beam and the one or moreauxiliary reflected beams and wherein the compounding function isdefined by the following equation:S _(out) =S _(main)−Σ_(m) W _(m) {tilde over (S)} _(aux) _(m) whereinS_(main) is the local measured signal for the main reflected beam;{tilde over (S)}_(aux) _(m) is the local measured signal for the m'thauxiliary beam, normalized so that the gain of all auxiliary receivebeams at an angular direction corresponding to the center of the mainreceive beam would match that of the main receive beam; S_(out) is theclutter suppressed local signal; and W_(m) is a local complex weight(clutter suppression weight) for the m'th auxiliary beam, which may becalculated for each voxel and/or each frame or a subset thereof.
 39. Themethod of claim 22, wherein the compounding function is computed by alinear combination of the outputs of another compounding function andthe main reflected beam signal, wherein the weights are based on anormalized generalized metric, and wherein the compounding function isdefined by the following equation:$S_{out} = {{\left\lbrack {1 - m} \right\rbrack S_{main}} + {m\frac{S_{main} - {W{\overset{\sim}{S}}_{aux}}}{1 - W}}}$wherein S_(main) is the local measured signal for the main reflectedbeam; {tilde over (S)}_(aux) is the local measured signal for anauxiliary beam, normalized so that the gain of the auxiliary beam at anangular direction corresponding to the center of the main receive beamwould match that of the main receive beam; m is a normalized generalizedmetric; S_(out) is the clutter suppressed local signal; and W is a localcomplex weight (clutter suppression weight), which may be calculated foreach voxel and/or each frame or a subset thereof.
 40. The method ofclaim 1, wherein the compounding function is computed by a linearcombination of the main reflected signal and the auxiliary reflectedsignal, defined by the following equation:$S_{out} = \frac{S_{main} - {W{\overset{\sim}{S}}_{aux}}}{1 - W}$wherein S_(main) is the local measured signal for the main reflectedbeam; {tilde over (S)}_(aux) is the local measured signal for anauxiliary beam, normalized so that the gain of the auxiliary beam at anangular direction corresponding to the center of the main receive beamwould match that of the main receive beam; and W is a local complexweight (clutter suppression weight), which may be calculated for eachvoxel and/or each frame or a subset thereof.
 41. The method of claim 39,wherein the clutter suppression weight for each applicable voxel {rightarrow over (p)} and/or each applicable frame n is estimated by thefollowing steps: (a) Define a region surrounding ({right arrow over(p)}, n), wherein the region may be spatial and/or temporal(spatial-temporal region); (b) Set the clutter suppression weight for({right arrow over (p)}, n) to be a function of at least one of: (i) thevalues of one or more generalized metrics and/or normalized generalizedmetrics within the spatial-temporal region; and (ii) the local complexratio between the signal for the main receive beam, and the signal inone of the normalized auxiliary receive beams, within thespatial-temporal region (local complex main-auxiliary ratio).
 42. Themethod of claim 41, wherein the clutter suppression weight for ({rightarrow over (p)}, n) is set to the local complex main-auxiliary ratio forthe voxel within the spatial-temporal region, whose generalized metricand/or normalized generalized metric values are most indicative ofsignificant clutter effects; wherein the computation of the cluttersuppression weight is based on the following assumptions: (a) Theclutter signal contribution changes slowly enough in space and/or intime, such that it is approximately constant within smallspatial-temporal regions; and (b) Speckle effects in ultrasonic imagingcause local variations in the clutter free signal levels even forapproximately homogeneous tissue regions, such that one may expect atleast one of the voxels in the spatial-temporal region to include lowfree signal levels.
 43. The method of claim 40, wherein a generalizedmetric is at least one of: (i) indicative of the probability for thecorresponding voxel to be substantially affected by clutter effectsonly, as derived using certain physical and/or physiologic assumptions;(ii) indicative of the probability for the corresponding voxel to besubstantially unaffected by clutter effects, as derived using certainphysical and/or physiologic assumptions; (iii) indicative of thepercentage of the received energy within the corresponding voxel thatoriginates from clutter effects; and (iv) set to a certain constant, ifthe corresponding voxel is not significantly affected by cluttereffects, and otherwise to a different constant, and wherein a normalizedgeneralized metric is a generalized metric defined as a real numberwhose value ranges from 0.0 to 1.0, and wherein the clutter suppressionweight for each applicable voxel and/or each applicable frame n isestimated by the following steps: (a) Define a region surrounding({right arrow over (p)}, n), wherein the region may be spatial and/ortemporal (spatial-temporal region); (b) Set the clutter suppressionweight for ({right arrow over (p)}, n) to be a function of at least oneof: (i) the values of one or more generalized metrics and/or normalizedgeneralized metrics within the spatial-temporal region; and (ii) thelocal complex ratio between the signal for the main receive beam, andthe signal in one of the normalized auxiliary receive beams, within thespatial-temporal region (local complex main-auxiliary ratio).
 44. Themethod of claim 43, wherein the clutter suppression weight for ({rightarrow over (p)}, n) is set to the local complex main-auxiliary ratio forthe voxel within the spatial-temporal region, whose generalized metricand/or normalized generalized metric values are most indicative ofsignificant clutter effects; wherein the computation of the cluttersuppression weight is based on the following assumptions: (a) Theclutter signal contribution changes slowly enough in space and/or intime, such that it is approximately constant within smallspatial-temporal regions; and (b) Speckle effects in ultrasonic imagingcause local variations in the clutter free signal levels even forapproximately homogeneous tissue regions, such that one may expect atleast one of the voxels in the spatial-temporal region to include lowfree signal levels.
 45. The method of claim 22, wherein two or morenormalized generalized metrics (m₁ and m₂) are defined and utilized bythe compounding function, defined by the following equation:S _(out)=[1−m ₁][1−m ₂ ]S _(main) +m ₁[1−m ₂ ]S ₁+[1−m ₁ ]m ₂ S ₂ +m ₁ m₂ S _(1,2)+ wherein S_(main) is the local measured signal for the mainreflected beam; S_(out) is the clutter suppressed local signal; S₁ andS₂ are the outputs of clutter suppression processes based on m₁ and m₂respectively; and S_(1,2) is the output of a clutter suppression processbased on both m₁ and m₂.
 46. The method of claim 22, wherein two or morenormalized generalized metrics are defined and utilized in thecompounding function, wherein at least one normalized generalized metricm₁ is used directly and at least one normalized generalized metric m₂ isused to compute a clutter suppression weight, wherein a cluttersuppression weight W₂ is computed using m₂, and the following equationmay be used:$S_{out} = {{{\left\lbrack {1 - m_{1}} \right\rbrack \left\lbrack {1 - m_{2}} \right\rbrack}S_{main}} + {m_{2}\frac{S_{main} - {W_{2}{\overset{\sim}{S}}_{aux}}}{1 - W_{2}}}}$wherein S_(main) is the local measured signal for the main reflectedbeam; and {tilde over (S)}_(aux) is the local measured signal for anauxiliary beam, normalized so that the gain of the auxiliary beam at anangular direction corresponding to the center of the main receive beamwould match that of the main receive beam.
 47. The method of claim 1,wherein the dataset collected by one or more reflected signals, in eachtime-swath or frame, is organized in a 1D, 2D or 3D scanned data array,wherein at least one of the following post-processing steps is appliedto the output of the compounding function: (a) For each applicable entryinto the scanned data array and/or each applicable frame, replace theoutput of the compounding function with the corresponding voxel of themain receive beam if a certain criterion is met; (b) For each applicableentry into the scanned data array and/or each applicable frame, replacethe local output of the compounding function with the mean signal,weighted mean signal or median signal within a small local region of thecompounding function output, defined in space and/or in time; and (c)For each applicable entry into the scanned data array and/or eachapplicable frame (current voxel), apply the following process: (i)Define a block within the scanned data array, approximately centered atthe current voxel (data block), wherein said block may beone-dimensional, two-dimensional, three dimensional or four-dimensional;(ii) For each of the main receive beam and the output of the compoundingfunction, calculate the mean magnitude or the weighted mean magnitudewithin the data block, including or excluding the current voxel, anddivide the result by the magnitude of the current voxel or by the meanmagnitude in a smaller block within the scanned data array; and (iii)Replace the value of the current voxel within the compounding functionoutput by one of: the value of the corresponding voxel in the mainreceive beam, the mean signal, weighted mean signal or median signalwithin the data block of the compounding function output; if the ratiobetween the output of step (ii) for the main receive beam and for theoutput of the compounding function is higher-than or lower-than than apredefined constant.
 48. The method of claim 1, wherein a plurality ofmain receive beams are used concurrently, and wherein at least one ofthe following: (a) A compounding model probability density function; (b)A generalized metric and/or a normalized generalized metric; (c) Aclutter suppression weight; and (d) The compounding function; iscomputed in at least one of the following ways: (i) Employing each pairof main receive beam and associated auxiliary receive beam used(main-auxiliary beam pair) separately; and (ii) Once for at least twomain-auxiliary beam pairs, wherein transformations are applied to datarelating to each main-auxiliary beam pair.
 49. The method of claim 1,wherein the input datasets from the main reflected signal and the one ormore auxiliary reflected signals are in one of the following forms: (a)Analog signal; or (b) Digital signal.
 50. The method of claim 1, whereinthe input datasets from the main reflected signal and the one or moreauxiliary reflected signals are in one of the following forms: (a) Realrepresentation; or (b) Complex representation.
 51. The method of claim1, wherein the input datasets from the main reflected signal and the oneor more auxiliary reflected signals are taken from at least one of thefollowing processing phases of ultrasound scanners: (a) Before matchedfiltering and before down-conversion; (b) Before matched filtering butafter down-conversion; (c) After matched filtering but beforedown-conversion; (d) After matched filtering and after down-conversion.